TFG JOANA

Carregar llibreries

library(openxlsx)
library(readxl)
library(dplyr)
library(vegan)
library(MASS)
library(reshape2)
library(reshape)
library(ggplot2)
library(eulerr)
library(purrr)
library(stringr)
library(ggbiplot)
library(indicspecies)
library(purrr)
library(lme4)
library(Matrix)
library(lmerTest)

Carregar la base de dades original (corregida) i eliminar les dades de Mordèl·lids i de les zones “Llig”, “Mura”, i “TV-7041”

# Obtenir el directori de l'arxiu actual
directorio <- getwd()

dades_finals <- read_excel(file.path(directorio,"Carreteres_ITS+altres.xlsx"))
dades_finals <- dades_finals %>%
  filter(Familia != "Mordellidae" & !Zona %in% c("Llig", "Mura", "TV-7041"))
head(dades_finals)
## # A tibble: 6 × 25
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6        28 Campus Juliol         7  2021      4 NS         <NA>         
## # ℹ 17 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>
# Guardar la nova base de dades
write.xlsx(dades_finals, file.path(directorio, "dades_finals.xlsx"))

ANÀLISI DESCRIPTIVA

# Definir nova variable: "Mostreig", segons Zona, Mes i Any
dades_finals <- dades_finals %>%
  mutate(Mostreig = paste(Zona, Mes, Any, sep = "-"))

Riquesa

Riquesa total

riquesa_total<- dades_finals %>% 
  summarise(riquesa_total=n_distinct(ID))
riquesa_total
## # A tibble: 1 × 1
##   riquesa_total
##           <int>
## 1           137

Riquesa per cada mostreig i tractament

riquesa <- dades_finals %>%
  group_by(Mostreig, Tractament) %>%
  summarise(riquesa = n_distinct(ID))
## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.
riquesa
## # A tibble: 11 × 3
## # Groups:   Mostreig [6]
##    Mostreig             Tractament riquesa
##    <chr>                <chr>        <int>
##  1 Campus-Juliol-2021   NS              20
##  2 Campus-Juliol-2021   S               21
##  3 Campus-Jun-2021      NS              45
##  4 Campus-Jun-2021      S               34
##  5 Campus-juliol-2020   NS              39
##  6 Campus-juliol-2020   S               24
##  7 Campus-maig-2020     NS              45
##  8 Franqueses-Maig-2021 NS              20
##  9 Franqueses-Maig-2021 S               11
## 10 Sabadell-Juny-2021   NS              16
## 11 Sabadell-Juny-2021   S               12

Riquesa per cada zona i tractament

riquesa_zones <- dades_finals %>%
  group_by(Zona, Tractament) %>%
  summarise(riquesa = n_distinct(ID))
## `summarise()` has grouped output by 'Zona'. You can override using the
## `.groups` argument.
riquesa_zones
## # A tibble: 6 × 3
## # Groups:   Zona [3]
##   Zona       Tractament riquesa
##   <chr>      <chr>        <int>
## 1 Campus     NS              98
## 2 Campus     S               59
## 3 Franqueses NS              20
## 4 Franqueses S               11
## 5 Sabadell   NS              16
## 6 Sabadell   S               12

Riquesa per grups

# Abelles
riquesa_abelles<- dades_finals %>% 
  filter(Ordre == 'HymenopteraAbella') %>% 
  summarise(riquesa_abelles=n_distinct(ID))
riquesa_abelles
## # A tibble: 1 × 1
##   riquesa_abelles
##             <int>
## 1              80
# Coleòpters
riquesa_coleoptera<- dades_finals %>% 
  filter(Ordre == 'Coleoptera') %>% 
  summarise(riquesa_coleoptera=n_distinct(ID))
riquesa_coleoptera
## # A tibble: 1 × 1
##   riquesa_coleoptera
##                <int>
## 1                 33
# Sírfids
riquesa_sirfids<- dades_finals %>% 
  filter(Familia == 'Syrphidae') %>% 
  summarise(riquesa_sirfids=n_distinct(ID))
riquesa_sirfids
## # A tibble: 1 × 1
##   riquesa_sirfids
##             <int>
## 1               9

Abundància

Abundància total

abundancia_total<-dades_finals %>% 
  summarise(Numero_total=n())
abundancia_total
## # A tibble: 1 × 1
##   Numero_total
##          <int>
## 1         1111

Abundància per cada zona i tractament

abundancia_zona <- dades_finals %>%
  group_by(Zona, Tractament) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Zona'. You can override using the
## `.groups` argument.
abundancia_zona
## # A tibble: 6 × 3
## # Groups:   Zona [3]
##   Zona       Tractament Numero_total
##   <chr>      <chr>             <int>
## 1 Campus     NS                  649
## 2 Campus     S                   323
## 3 Franqueses NS                   45
## 4 Franqueses S                    14
## 5 Sabadell   NS                   29
## 6 Sabadell   S                    51

Abundància de cada gènere total

ab_genere <- dades_finals %>%
  group_by(Genere) %>%
  summarise(Numero_total = n())
head(ab_genere)
## # A tibble: 6 × 2
##   Genere     Numero_total
##   <chr>             <int>
## 1 Acmaeodera            6
## 2 Amegilla              1
## 3 Andrena               6
## 4 Anthaxia             29
## 5 Anthidium            10
## 6 Apis                 11

Abundància de cada ID per zona i tractament

abundaciazona_ID <- dades_finals %>%
  group_by(Zona, Tractament, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Zona', 'Tractament'. You can override
## using the `.groups` argument.
head(abundaciazona_ID)
## # A tibble: 6 × 4
## # Groups:   Zona, Tractament [1]
##   Zona   Tractament ID                     Numero_total
##   <chr>  <chr>      <chr>                         <int>
## 1 Campus NS         Acmaeodera cylindrica             5
## 2 Campus NS         Amegilla sp.                      1
## 3 Campus NS         Anthaxia (Anthaxia)               1
## 4 Campus NS         Anthaxia (Melanthaxia)            5
## 5 Campus NS         Anthaxia hungarica               10
## 6 Campus NS         Anthaxia sp.                      7

Abundancia de cada ID per cada mostreig i tractament

abundacia_ID <- dades_finals %>%
  group_by(Mostreig, Tractament, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Mostreig', 'Tractament'. You can override
## using the `.groups` argument.
head(abundacia_ID)
## # A tibble: 6 × 4
## # Groups:   Mostreig, Tractament [1]
##   Mostreig           Tractament ID                   Numero_total
##   <chr>              <chr>      <chr>                       <int>
## 1 Campus-Juliol-2021 NS         Amegilla sp.                    1
## 2 Campus-Juliol-2021 NS         Anthaxia sp.                    1
## 3 Campus-Juliol-2021 NS         Bombus pascuorum                1
## 4 Campus-Juliol-2021 NS         Ceratina cucurbitina            1
## 5 Campus-Juliol-2021 NS         Ceratina parvula                7
## 6 Campus-Juliol-2021 NS         Chrysididae                     1

DIVERSITAT DE SHANNON I CORBES D’ACUMULACIÓ D’ESPÈCIES

Diversitat de Shannon per Zona i Tractament

# Preparar dades: Crear una taula de l'abundància de cada ID (columnes) segons la zona i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux2 <- cast(abundaciazona_ID, Zona + Tractament ~ ID, value='Numero_total', FUN=mean)
taula_aux2 <- as.data.frame(taula_aux2)
taula_aux2[is.na(taula_aux2)] <- 0
# Agrupar Zona amb Tractament (en una nova variable Agrupacio2)
data_shannon_ZT <- taula_aux2 %>%
  mutate(Agrupacio2 = paste(Zona, Tractament, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Zona, -Tractament)
# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio2 i establir-la com a identificador de les files
rownames(data_shannon_ZT) <- data_shannon_ZT$Agrupacio2
data_shannon_ZT <- data_shannon_ZT %>% dplyr::select(-Agrupacio2) %>% mutate_all(as.numeric)
# Càlcul de la diversitat de Shannon per cada Mostreig i Tractament
shannondiv2 <- diversity(data_shannon_ZT)
print(shannondiv2)
##     Campus_NS      Campus_S Franqueses_NS  Franqueses_S   Sabadell_NS 
##      3.612943      2.870658      2.627033      2.341994      2.382088 
##    Sabadell_S 
##      2.003690

Diversitat de Shannon per Mostreig i Tractament

# Preparar dades: Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux <- cast(abundacia_ID, Mostreig + Tractament ~ ID, value='Numero_total', FUN=mean)
taula_aux <- as.data.frame(taula_aux)
taula_aux[is.na(taula_aux)] <- 0
# Agrupar Mostreig amb Tractament (en una nova variable Agrupacio)
data_shannon_MT <- taula_aux %>%
  mutate(Agrupacio = paste(Mostreig, Tractament, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament)
# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_MT) <- data_shannon_MT$Agrupacio
data_shannon_MT <- data_shannon_MT %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)
# Càlcul de la diversitat de Shannon per cada Mostreig i Tractament
shannondiv <- diversity(data_shannon_MT)
print(shannondiv)
##   Campus-juliol-2020_NS    Campus-juliol-2020_S   Campus-Juliol-2021_NS 
##                3.011325                2.404801                2.593731 
##    Campus-Juliol-2021_S      Campus-Jun-2021_NS       Campus-Jun-2021_S 
##                1.933200                3.161179                2.856531 
##     Campus-maig-2020_NS Franqueses-Maig-2021_NS  Franqueses-Maig-2021_S 
##                2.957046                2.627033                2.341994 
##   Sabadell-Juny-2021_NS    Sabadell-Juny-2021_S 
##                2.382088                2.003690

Corba d’acumulació d’espècies per mostreig (11 mostrejos)

corba <- specaccum(data_shannon_MT, permutations = 500)
## Warning in cor(x > 0): the standard deviation is zero
plot(corba)

Extrapolació de dades (quantes espècies més trobaríem si féssim més mostrejos)

specpool(data_shannon_MT)
##     Species     chao  chao.se jack1 jack1.se    jack2     boot  boot.se  n
## All     137 236.8148 31.40658   207 26.07681 250.1182 167.3704 12.88879 11
plot(poolaccum(data_shannon_MT))

En aquest mostreig s’han trobat 137 espècies (riquesa), però s’estima que hi ha altres espècies “unseen” (no detectades). Segons diferents estimacions, la riquesa total tenint en compte aquestes espècies no detectades seria 236 (chao), 207 (jackknife primer ordre), 250 (jackknife segon ordre), i 167 (bootstrap).

Diversitat de Shannon per Mostreig, Tractament i Trampa

# Preparar  dades:
#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa
abundacia_trampa <- dades_finals %>%
  group_by(Mostreig, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.
abundacia_trampa
## # A tibble: 564 × 5
## # Groups:   Mostreig, Tractament, Trampa [101]
##    Mostreig           Tractament Trampa ID                         Numero_total
##    <chr>              <chr>       <dbl> <chr>                             <int>
##  1 Campus-Juliol-2021 NS              2 Lasioglossum glabriusculum            1
##  2 Campus-Juliol-2021 NS              4 Curculionidae                         1
##  3 Campus-Juliol-2021 NS              4 Lasioglossum malachurum               1
##  4 Campus-Juliol-2021 NS              4 Pompilidae                            1
##  5 Campus-Juliol-2021 NS              4 Pyronia cecilia                       1
##  6 Campus-Juliol-2021 NS              6 Lasioglossum glabriusculum            1
##  7 Campus-Juliol-2021 NS              6 Lasioglossum politum                  3
##  8 Campus-Juliol-2021 NS              6 Pompilidae                            1
##  9 Campus-Juliol-2021 NS              8 Chrysididae                           1
## 10 Campus-Juliol-2021 NS              8 Hylaeus sp.                           1
## # ℹ 554 more rows
#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux <- cast(abundacia_trampa, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_aux <- as.data.frame(taula_aux)
taula_aux[is.na(taula_aux)] <- 0
# Agrupar Mostreig amb Tractament i Trampa
data_shannon_MTTR <- taula_aux %>%
  mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament, -Trampa)
# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_MTTR) <- data_shannon_MTTR$Agrupacio
data_shannon_MTTR <- data_shannon_MTTR %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)
# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondiv <- diversity(data_shannon_MTTR)
print(shannondiv)
##   Campus-juliol-2020_NS_2   Campus-juliol-2020_NS_3   Campus-juliol-2020_NS_4 
##                 1.7917595                 1.8310205                 2.0449312 
##   Campus-juliol-2020_NS_6   Campus-juliol-2020_NS_8  Campus-juliol-2020_NS_10 
##                 1.7328680                 1.8310205                 1.0735428 
##  Campus-juliol-2020_NS_11  Campus-juliol-2020_NS_12  Campus-juliol-2020_NS_13 
##                 0.0000000                 1.7478681                 1.6434177 
##  Campus-juliol-2020_NS_14  Campus-juliol-2020_NS_15  Campus-juliol-2020_NS_16 
##                 1.1240357                 1.5498260                 1.3862944 
##  Campus-juliol-2020_NS_20    Campus-juliol-2020_S_1    Campus-juliol-2020_S_5 
##                 2.2718685                 0.6931472                 1.9072840 
##    Campus-juliol-2020_S_7    Campus-juliol-2020_S_9   Campus-juliol-2020_S_17 
##                 0.6924212                 0.0000000                 1.0986123 
##   Campus-juliol-2020_S_18   Campus-juliol-2020_S_19   Campus-juliol-2020_S_21 
##                 1.6094379                 1.5607104                 1.1189689 
##   Campus-juliol-2020_S_22   Campus-Juliol-2021_NS_2   Campus-Juliol-2021_NS_4 
##                 1.0751393                 0.0000000                 1.3862944 
##   Campus-Juliol-2021_NS_6   Campus-Juliol-2021_NS_8  Campus-Juliol-2021_NS_10 
##                 0.9502705                 1.4648164                 0.6931472 
##  Campus-Juliol-2021_NS_12  Campus-Juliol-2021_NS_14  Campus-Juliol-2021_NS_16 
##                 1.3296613                 1.3321790                 0.0000000 
##  Campus-Juliol-2021_NS_18  Campus-Juliol-2021_NS_20    Campus-Juliol-2021_S_1 
##                 1.3862944                 1.0986123                 1.0986123 
##    Campus-Juliol-2021_S_3    Campus-Juliol-2021_S_5    Campus-Juliol-2021_S_7 
##                 0.6365142                 2.0963526                 1.5000764 
##    Campus-Juliol-2021_S_9   Campus-Juliol-2021_S_11   Campus-Juliol-2021_S_13 
##                 1.1973401                 1.4388830                 0.0000000 
##   Campus-Juliol-2021_S_15   Campus-Juliol-2021_S_17      Campus-Jun-2021_NS_2 
##                 0.6931472                 0.5091373                 2.3345491 
##      Campus-Jun-2021_NS_4      Campus-Jun-2021_NS_6      Campus-Jun-2021_NS_8 
##                 2.0854684                 1.1537419                 2.6197294 
##     Campus-Jun-2021_NS_10     Campus-Jun-2021_NS_12     Campus-Jun-2021_NS_14 
##                 1.8310205                 2.5575077                 1.7623137 
##     Campus-Jun-2021_NS_16     Campus-Jun-2021_NS_18     Campus-Jun-2021_NS_20 
##                 1.3593685                 1.9061547                 0.9502705 
##       Campus-Jun-2021_S_1       Campus-Jun-2021_S_3       Campus-Jun-2021_S_5 
##                 2.2718685                 2.0028830                 1.7351265 
##       Campus-Jun-2021_S_7       Campus-Jun-2021_S_9      Campus-Jun-2021_S_11 
##                 0.9502705                 2.1383972                 1.5595812 
##      Campus-Jun-2021_S_13      Campus-Jun-2021_S_15      Campus-Jun-2021_S_17 
##                 1.8576598                 0.7356219                 2.1639557 
##      Campus-Jun-2021_S_19     Campus-maig-2020_NS_1     Campus-maig-2020_NS_2 
##                 1.9722470                 0.8392696                 1.3642812 
##     Campus-maig-2020_NS_3     Campus-maig-2020_NS_4     Campus-maig-2020_NS_5 
##                 1.6957425                 2.5521719                 1.5171064 
##     Campus-maig-2020_NS_6     Campus-maig-2020_NS_7     Campus-maig-2020_NS_8 
##                 2.2614951                 1.2424533                 1.7480673 
##     Campus-maig-2020_NS_9    Campus-maig-2020_NS_10    Campus-maig-2020_NS_11 
##                 1.6769878                 1.9722470                 1.7480673 
##    Campus-maig-2020_NS_12    Campus-maig-2020_NS_13    Campus-maig-2020_NS_14 
##                 1.0986123                 1.5481155                 1.8891592 
##    Campus-maig-2020_NS_15    Campus-maig-2020_NS_16    Campus-maig-2020_NS_17 
##                 1.8392967                 1.0356097                 1.3296613 
##    Campus-maig-2020_NS_18    Campus-maig-2020_NS_19    Campus-maig-2020_NS_20 
##                 1.5810938                 1.3535914                 0.9002561 
## Franqueses-Maig-2021_NS_1 Franqueses-Maig-2021_NS_2 Franqueses-Maig-2021_NS_3 
##                 1.0789922                 0.8675632                 0.6931472 
## Franqueses-Maig-2021_NS_4 Franqueses-Maig-2021_NS_5  Franqueses-Maig-2021_S_1 
##                 2.2739657                 0.7549968                 1.3321790 
##  Franqueses-Maig-2021_S_2  Franqueses-Maig-2021_S_4   Sabadell-Juny-2021_NS_1 
##                 1.0986123                 1.5607104                 1.7478681 
##   Sabadell-Juny-2021_NS_4   Sabadell-Juny-2021_NS_7   Sabadell-Juny-2021_NS_8 
##                 0.9502705                 1.7917595                 1.4750763 
##  Sabadell-Juny-2021_NS_11    Sabadell-Juny-2021_S_2    Sabadell-Juny-2021_S_3 
##                 1.3862944                 0.6365142                 0.5623351 
##    Sabadell-Juny-2021_S_5    Sabadell-Juny-2021_S_6    Sabadell-Juny-2021_S_9 
##                 1.4995094                 1.0986123                 0.6931472 
##   Sabadell-Juny-2021_S_10   Sabadell-Juny-2021_S_12 
##                 0.5004024                 1.2852930

Corba d’acumulació d’espècies segons mostreig, tractament i trampa

corba <- specaccum(data_shannon_MTTR, permutations = 500)
plot(corba)

Extrapolació de dades (quantes espècies més trobaríem si féssim més mostrejos)

specpool(data_shannon_MTTR)
##     Species     chao  chao.se    jack1 jack1.se    jack2     boot  boot.se   n
## All     137 224.1493 29.51259 201.3564 11.27322 241.7798 164.5183 5.987207 101
plot(poolaccum(data_shannon_MTTR))

Amb aquest mostreig s’han trobat 137 espècies (riquesa), però s’estima que hi ha altres espècies “unseen” (no detectades). Segons diferents estimacions, la riquesa total tenint en compte aquestes espècies no detectades seria 224 (chao), 201 (jackknife primer ordre), 241 (jackknife segon ordre), i 164 (bootstrap).

Corba i extrapolació pel Campus:

#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa, amb dades filtrades pel Campus
abundacia_trampa_campus <- dades_finals %>%
  filter(Zona == 'Campus') %>%
  group_by(Mostreig, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.
#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_auxcampus <- cast(abundacia_trampa_campus, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_auxcampus <- as.data.frame(taula_auxcampus)
taula_auxcampus[is.na(taula_auxcampus)] <- 0

# Agrupar Mostreig amb Tractament i Trampa
data_shannon_campus <- taula_auxcampus %>%
  mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament, -Trampa)

# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_campus) <- data_shannon_campus$Agrupacio
data_shannon_campus <- data_shannon_campus %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)

# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondivcampus <- diversity(data_shannon_campus)
print(shannondivcampus)
##  Campus-juliol-2020_NS_2  Campus-juliol-2020_NS_3  Campus-juliol-2020_NS_4 
##                1.7917595                1.8310205                2.0449312 
##  Campus-juliol-2020_NS_6  Campus-juliol-2020_NS_8 Campus-juliol-2020_NS_10 
##                1.7328680                1.8310205                1.0735428 
## Campus-juliol-2020_NS_11 Campus-juliol-2020_NS_12 Campus-juliol-2020_NS_13 
##                0.0000000                1.7478681                1.6434177 
## Campus-juliol-2020_NS_14 Campus-juliol-2020_NS_15 Campus-juliol-2020_NS_16 
##                1.1240357                1.5498260                1.3862944 
## Campus-juliol-2020_NS_20   Campus-juliol-2020_S_1   Campus-juliol-2020_S_5 
##                2.2718685                0.6931472                1.9072840 
##   Campus-juliol-2020_S_7   Campus-juliol-2020_S_9  Campus-juliol-2020_S_17 
##                0.6924212                0.0000000                1.0986123 
##  Campus-juliol-2020_S_18  Campus-juliol-2020_S_19  Campus-juliol-2020_S_21 
##                1.6094379                1.5607104                1.1189689 
##  Campus-juliol-2020_S_22  Campus-Juliol-2021_NS_2  Campus-Juliol-2021_NS_4 
##                1.0751393                0.0000000                1.3862944 
##  Campus-Juliol-2021_NS_6  Campus-Juliol-2021_NS_8 Campus-Juliol-2021_NS_10 
##                0.9502705                1.4648164                0.6931472 
## Campus-Juliol-2021_NS_12 Campus-Juliol-2021_NS_14 Campus-Juliol-2021_NS_16 
##                1.3296613                1.3321790                0.0000000 
## Campus-Juliol-2021_NS_18 Campus-Juliol-2021_NS_20   Campus-Juliol-2021_S_1 
##                1.3862944                1.0986123                1.0986123 
##   Campus-Juliol-2021_S_3   Campus-Juliol-2021_S_5   Campus-Juliol-2021_S_7 
##                0.6365142                2.0963526                1.5000764 
##   Campus-Juliol-2021_S_9  Campus-Juliol-2021_S_11  Campus-Juliol-2021_S_13 
##                1.1973401                1.4388830                0.0000000 
##  Campus-Juliol-2021_S_15  Campus-Juliol-2021_S_17     Campus-Jun-2021_NS_2 
##                0.6931472                0.5091373                2.3345491 
##     Campus-Jun-2021_NS_4     Campus-Jun-2021_NS_6     Campus-Jun-2021_NS_8 
##                2.0854684                1.1537419                2.6197294 
##    Campus-Jun-2021_NS_10    Campus-Jun-2021_NS_12    Campus-Jun-2021_NS_14 
##                1.8310205                2.5575077                1.7623137 
##    Campus-Jun-2021_NS_16    Campus-Jun-2021_NS_18    Campus-Jun-2021_NS_20 
##                1.3593685                1.9061547                0.9502705 
##      Campus-Jun-2021_S_1      Campus-Jun-2021_S_3      Campus-Jun-2021_S_5 
##                2.2718685                2.0028830                1.7351265 
##      Campus-Jun-2021_S_7      Campus-Jun-2021_S_9     Campus-Jun-2021_S_11 
##                0.9502705                2.1383972                1.5595812 
##     Campus-Jun-2021_S_13     Campus-Jun-2021_S_15     Campus-Jun-2021_S_17 
##                1.8576598                0.7356219                2.1639557 
##     Campus-Jun-2021_S_19    Campus-maig-2020_NS_1    Campus-maig-2020_NS_2 
##                1.9722470                0.8392696                1.3642812 
##    Campus-maig-2020_NS_3    Campus-maig-2020_NS_4    Campus-maig-2020_NS_5 
##                1.6957425                2.5521719                1.5171064 
##    Campus-maig-2020_NS_6    Campus-maig-2020_NS_7    Campus-maig-2020_NS_8 
##                2.2614951                1.2424533                1.7480673 
##    Campus-maig-2020_NS_9   Campus-maig-2020_NS_10   Campus-maig-2020_NS_11 
##                1.6769878                1.9722470                1.7480673 
##   Campus-maig-2020_NS_12   Campus-maig-2020_NS_13   Campus-maig-2020_NS_14 
##                1.0986123                1.5481155                1.8891592 
##   Campus-maig-2020_NS_15   Campus-maig-2020_NS_16   Campus-maig-2020_NS_17 
##                1.8392967                1.0356097                1.3296613 
##   Campus-maig-2020_NS_18   Campus-maig-2020_NS_19   Campus-maig-2020_NS_20 
##                1.5810938                1.3535914                0.9002561
corbacampus <- specaccum(data_shannon_campus, permutations = 500)
plot(corbacampus)

# Extrapolació de les dades
specpool(data_shannon_campus)
##     Species     chao  chao.se    jack1 jack1.se    jack2     boot  boot.se  n
## All     118 184.0553 24.55884 170.3457 9.341038 201.8116 140.7347 5.098803 81

Corba i extrapolació per Franqueses:

#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa, amb dades filtrades per Franqueses
abundacia_trampa_franqueses <- dades_finals %>%
  filter(Zona == 'Franqueses') %>%
  group_by(Mostreig, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.
#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_auxfranq <- cast(abundacia_trampa_franqueses, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_auxfranq <- as.data.frame(taula_auxfranq)
taula_auxfranq[is.na(taula_auxfranq)] <- 0

# Agrupar Mostreig amb Tractament i Trampa
data_shannon_franqueses <- taula_auxfranq %>%
  mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament, -Trampa)

# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_franqueses) <- data_shannon_franqueses$Agrupacio
data_shannon_franqueses <- data_shannon_franqueses %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)

# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondivfranq <- diversity(data_shannon_franqueses)
print(shannondivfranq)
## Franqueses-Maig-2021_NS_1 Franqueses-Maig-2021_NS_2 Franqueses-Maig-2021_NS_3 
##                 1.0789922                 0.8675632                 0.6931472 
## Franqueses-Maig-2021_NS_4 Franqueses-Maig-2021_NS_5  Franqueses-Maig-2021_S_1 
##                 2.2739657                 0.7549968                 1.3321790 
##  Franqueses-Maig-2021_S_2  Franqueses-Maig-2021_S_4 
##                 1.0986123                 1.5607104
corbafranqueses <- specaccum(data_shannon_franqueses, permutations = 500)
plot(corbafranqueses)

# Extrapolació de les dades
specpool(data_shannon_franqueses)
##     Species     chao  chao.se  jack1 jack1.se    jack2     boot  boot.se n
## All      26 52.32292 16.96108 42.625 8.605049 53.01786 33.13315 4.330973 8

Corba i extrapolació per Sabadell

#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa, amb dades filtrades per Sabadell
abundacia_trampa_sabadell <- dades_finals %>%
  filter(Zona == 'Sabadell') %>%
  group_by(Mostreig, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.
#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_auxsaba <- cast(abundacia_trampa_sabadell, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_auxsaba <- as.data.frame(taula_auxsaba)
taula_auxsaba[is.na(taula_auxsaba)] <- 0

# Agrupar Mostreig amb Tractament i Trampa
data_shannon_sabadell <- taula_auxsaba %>%
  mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament, -Trampa)

# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_sabadell) <- data_shannon_sabadell$Agrupacio
data_shannon_sabadell <- data_shannon_sabadell %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)

# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondivsaba <- diversity(data_shannon_sabadell)
print(shannondivsaba)
##  Sabadell-Juny-2021_NS_1  Sabadell-Juny-2021_NS_4  Sabadell-Juny-2021_NS_7 
##                1.7478681                0.9502705                1.7917595 
##  Sabadell-Juny-2021_NS_8 Sabadell-Juny-2021_NS_11   Sabadell-Juny-2021_S_2 
##                1.4750763                1.3862944                0.6365142 
##   Sabadell-Juny-2021_S_3   Sabadell-Juny-2021_S_5   Sabadell-Juny-2021_S_6 
##                0.5623351                1.4995094                1.0986123 
##   Sabadell-Juny-2021_S_9  Sabadell-Juny-2021_S_10  Sabadell-Juny-2021_S_12 
##                0.6931472                0.5004024                1.2852930
corbasabadell <- specaccum(data_shannon_sabadell, permutations = 500)
plot(corbasabadell)

# Extrapolació de les dades
specpool(data_shannon_sabadell)
##     Species     chao  chao.se    jack1 jack1.se   jack2     boot  boot.se  n
## All      23 45.45833 17.10731 35.83333 5.316379 44.4697 28.44918 2.631169 12

T-student per comparar la diversitat de Shannon entre mostrejos i tractaments

#Convertir la taula de shannondiv en dataframe
shannon_ttest <- as.data.frame(shannondiv)

# Convertir l'ID de les files (mostreig + tractament + número de trampa) en una variable (nova columna)
shannon_ttest$row <- rownames(shannon_ttest)

# De la nova variable, separar en dues columnes el mostreig i eltractament (i no utilitzem el número de trampa)
shannon_ttest <- shannon_ttest %>% mutate(Tractament = str_split(row, "_", simplify = TRUE)[, 2], Mostreig = str_split(row, "_", simplify = TRUE)[,1]) %>% dplyr::select(-row)

# Prova t-student i boxplots
shannonT <- shannon_ttest %>% 
  group_by(Mostreig) %>% 
  summarise(
    mean_S = mean(shannondiv[Tractament=='S']),
    mean_NS = mean(shannondiv[Tractament=='NS']),
    t_test_result = ifelse(length(unique(Tractament)) == 2, 
                             t.test(shannondiv ~ Tractament)$p.value,
                             NA),.groups = "drop")
  
shannonT
## # A tibble: 6 × 4
##   Mostreig              mean_S mean_NS t_test_result
##   <chr>                  <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021     1.02    0.964        0.844 
## 2 Campus-Jun-2021        1.74    1.86         0.637 
## 3 Campus-juliol-2020     1.08    1.54         0.0840
## 4 Campus-maig-2020     NaN       1.56        NA     
## 5 Franqueses-Maig-2021   1.33    1.13         0.566 
## 6 Sabadell-Juny-2021     0.897   1.47         0.0232
ggplot(shannon_ttest, aes(x = factor(Mostreig), y = shannondiv, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Diversitat de Shannon x Zona i Tractaments",
     x = "Mostreig",
     y = "Diversitat Shannon")

Només en el mostreig de Sabadell la diversitat de Shannon és significativament diferent (p valor < 0.05): l’índex de Shannon és major al tractament NS.

Diversitat de Shannon d’abelles (per mostreig)

# Preparar les dades: crear taula de dades d'abundància només per abelles
abundanciabelles_ID <- dades_finals %>% 
  mutate(Mostreig = paste(Zona, Mes, Any, sep = "-")) %>% 
  filter(Ordre == 'HymenopteraAbella') %>% 
  group_by(Mostreig, Tractament, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Mostreig', 'Tractament'. You can override
## using the `.groups` argument.
# View(abundanciabelles_ID)

#Crear una taula auxiliar de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux <- cast(abundanciabelles_ID, Mostreig + Tractament ~ ID, value='Numero_total', FUN=mean)
taula_aux <- as.data.frame(taula_aux)
taula_aux[is.na(taula_aux)] <- 0

#Agrupar Mostreig i tractament en una nova variable (Agrupacio)
data_shannon_abelles <- taula_aux %>% 
  mutate(Agrupacio = paste(Mostreig, Tractament, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament)

#Definir els rownames amb la variable Agrupacio (mostreig i tractament)
rownames(data_shannon_abelles) <- data_shannon_abelles$Agrupacio
data_shannon_abelles <- data_shannon_abelles %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)

#View(data_shannon_abelles)
# Càlcul diversitat de Shannon
shannondiv_abelles <- diversity(data_shannon_abelles)
print(shannondiv_abelles)
##   Campus-juliol-2020_NS    Campus-juliol-2020_S   Campus-Juliol-2021_NS 
##                2.611229                1.905899                2.090031 
##    Campus-Juliol-2021_S      Campus-Jun-2021_NS       Campus-Jun-2021_S 
##                1.806531                2.701045                2.080060 
##     Campus-maig-2020_NS Franqueses-Maig-2021_NS  Franqueses-Maig-2021_S 
##                2.199104                2.239746                1.560710 
##   Sabadell-Juny-2021_NS    Sabadell-Juny-2021_S 
##                2.045357                1.884985

Corbes d’acumulació d’espècies per mostreig

Corba d’acumulació. Riquesa en funció del nº de trampes (specaccum)

# A la taula data_shannon_MTTR extreure una nova columna que indiqui el mostreig
 metadata_shannon <- data_shannon_MTTR %>%
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) 

# View(metadata_shannon)

# Subset each mostreig into its own dataframe (i eliminar la columna de "mostreig" perquè totes les dades de la taula siguin numèriques)
metadata_shannon %>% filter(Zona == 'Campus') %>% dplyr::select(-Zona)-> Campus
metadata_shannon %>% filter(Zona == 'Franqueses') %>% dplyr::select(-Zona) -> Franqueses
metadata_shannon %>% filter(Zona == 'Llig') %>% dplyr::select(-Zona) -> Llig
metadata_shannon %>% filter(Zona == 'Mura') %>% dplyr::select(-Zona) -> Mura
metadata_shannon %>% filter(Zona == 'Sabadell') %>% dplyr::select(-Zona) -> Sabadell
metadata_shannon %>% filter(Zona == 'TV') %>% dplyr::select(-Zona) -> TV7041

# Corba d'acumulació d'espècies per cada mostreig
corba1 = specaccum(Campus, method = "random")
corba2 = specaccum(Franqueses, method = "random")
corba3 = specaccum(Sabadell, method = "random")

# Creating data frames for each set of data
data1 <- data.frame(Sites = corba1$sites, Richness = corba1$richness, SD = corba1$sd)
data1$Zona <- "Campus"

data2 <- data.frame(Sites = corba2$sites, Richness = corba2$richness, SD = corba2$sd)
data2$Zona <- "Franqueses"

data3 <- data.frame(Sites = corba3$sites, Richness = corba3$richness, SD = corba3$sd)
data3$Zona <- "Sabadell"


# Combine all data frames into one
dades_combinades <- rbind(data1, data2, data3)

library(ggplot2)

# Gráfico de dispersión con líneas y error bars
ggplot(dades_combinades, aes(x = Sites, y = Richness, color = Zona)) +
  geom_point() +
  geom_line(aes(group = Zona), alpha = 0.7) +  # Líneas
  # geom_ribbon(aes(ymin = Richness - 2 * SD, ymax = Richness + 2 * SD, fill = Mostreig), alpha = 0.01) +  # Incertidumbre
  geom_errorbar(aes(ymin = Richness - 2 * SD, ymax = Richness + 2 * SD), width = 0.2, alpha = 0.5) +  # Barras de error
  ylim(0,NA)+
  labs(x = "Nº trampes", y = "Riquesa", color = "Zona") +
  theme_minimal()

Corba de rarefracció. Riquesa en funció del nombre d’individus (rarefy)

dades_agrupades <- metadata_shannon %>%
  group_by(Zona) %>%
  summarise(across(everything(), sum)) %>% 
  as.data.frame()

#View (dades_agrupades) 

# Buscar el nombre mínim d'observacions en una zona (és 84, a TV juliol)
abundacia_trampa <- dades_finals %>%
  group_by(Zona, Mes, Any, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Zona', 'Mes', 'Any', 'Tractament',
## 'Trampa'. You can override using the `.groups` argument.
abundacia_trampa
## # A tibble: 564 × 7
## # Groups:   Zona, Mes, Any, Tractament, Trampa [101]
##    Zona   Mes      Any Tractament Trampa ID                         Numero_total
##    <chr>  <chr>  <dbl> <chr>       <dbl> <chr>                             <int>
##  1 Campus Juliol  2021 NS              2 Lasioglossum glabriusculum            1
##  2 Campus Juliol  2021 NS              4 Curculionidae                         1
##  3 Campus Juliol  2021 NS              4 Lasioglossum malachurum               1
##  4 Campus Juliol  2021 NS              4 Pompilidae                            1
##  5 Campus Juliol  2021 NS              4 Pyronia cecilia                       1
##  6 Campus Juliol  2021 NS              6 Lasioglossum glabriusculum            1
##  7 Campus Juliol  2021 NS              6 Lasioglossum politum                  3
##  8 Campus Juliol  2021 NS              6 Pompilidae                            1
##  9 Campus Juliol  2021 NS              8 Chrysididae                           1
## 10 Campus Juliol  2021 NS              8 Hylaeus sp.                           1
## # ℹ 554 more rows
min_n <- abundacia_trampa %>% 
  group_by(Zona) %>% 
  summarize (n_especimens = sum(Numero_total)) %>% 
  summarize(min=min(n_especimens)) %>% 
  pull(min)

rownames(dades_agrupades) <- dades_agrupades$Zona
dades_agrupades <- dades_agrupades[,-1]
dades_agrupades <- dades_agrupades %>% 
                    mutate_all(as.numeric)

vegan::rarefy(dades_agrupades,sample=min_n)
##     Campus Franqueses   Sabadell 
##   28.12144   26.00000   19.50913 
## attr(,"Subsample")
## [1] 59
rarecurve(dades_agrupades,step=5)

# Mirar vídeo i posar-ho més maco?

Corba d’interpolació: riquesa en funció del nombre d’individus (iNEXT). (??)

library(iNEXT)
transposades <- t(dades_agrupades)
D_abund <- iNEXT (transposades, datatype = 'abundance')
plot (D_abund)

Corbes de rarefracció d’abelles per cada tractament

# Extreure la columna de tractament a la taula d'abundàncies per ID de cada mostreig
dades_shannon_abelles <- data_shannon_abelles %>% 
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2])

# Agrupar per tractament
dades_agrupades_tractament <- dades_shannon_abelles %>%
  group_by(Tractament) %>%
  summarise(across(everything(), sum)) %>% 
  as.data.frame()

#View (dades_agrupades_tractament) 

dades_abelles <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]

# Buscar el nombre mínim d'observacions en una zona
abundacia_trampa_abelles <- dades_abelles %>%
  group_by(Zona, Mes, Any, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Zona', 'Mes', 'Any', 'Tractament',
## 'Trampa'. You can override using the `.groups` argument.
#View (abundacia_trampa_abelles)

min_n <- abundacia_trampa_abelles %>% 
  group_by(Tractament) %>% 
  summarize (n_especimens = sum(Numero_total)) %>% 
  summarize(min=min(n_especimens)) %>% 
  pull(min)

rownames(dades_agrupades_tractament) <- dades_agrupades_tractament$Tractament
dades_agrupades_tractament <- dades_agrupades_tractament[,-1]
dades_agrupades_tractament <- dades_agrupades_tractament %>% 
                    mutate_all(as.numeric)

vegan::rarefy(dades_agrupades_tractament,sample=min_n)
##       NS        S 
## 59.89088 37.00000 
## attr(,"Subsample")
## [1] 321
rarecurve(dades_agrupades_tractament,step=5)

Corbes de rarefracció d’abelles al Campus per tractament

# Extreure la columna Zona de dades_shannon_abelles i filtrar per Campus
dades_shannon_abellescampus <- dades_shannon_abelles %>% 
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1])
dades_shannon_abellescampus <- dades_shannon_abellescampus[dades_shannon_abellescampus$Zona == "Campus", ]

# Agrupar per tractament
dadescampus_agrupades_tractament <- dades_shannon_abellescampus %>%
  group_by(Tractament) %>%
  summarise(across(where(is.numeric), sum)) %>% 
  as.data.frame()

#View (dadescampus_agrupades_tractament) 

# Buscar el nombre mínim d'observacions en una zona
abundacia_trampa_abellescampus <- dades_abelles %>%
  group_by(Zona, Mes, Any, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Zona', 'Mes', 'Any', 'Tractament',
## 'Trampa'. You can override using the `.groups` argument.
#View (abundacia_trampa_abellescampus)

min_n <- abundacia_trampa_abellescampus %>% 
  group_by(Tractament) %>% 
  summarize (n_especimens = sum(Numero_total)) %>% 
  summarize(min=min(n_especimens)) %>% 
  pull(min)

rownames(dadescampus_agrupades_tractament) <- dadescampus_agrupades_tractament$Tractament
dadescampus_agrupades_tractament <- dadescampus_agrupades_tractament[,-1]
dadescampus_agrupades_tractament <- dadescampus_agrupades_tractament %>% 
                    mutate_all(as.numeric)

vegan::rarefy(dadescampus_agrupades_tractament,sample=min_n)
## Warning in vegan::rarefy(dadescampus_agrupades_tractament, sample = min_n):
## requested 'sample' was larger than smallest site maximum (266)
##       NS        S 
## 54.01949 35.00000 
## attr(,"Subsample")
## [1] 321
rarecurve(dadescampus_agrupades_tractament,step=5)

PROVES T-STUDENT

Comparació de la riquesa per mostrejos i tractaments

Riquesa total

riquesa_x_trampa <- dades_finals %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Riquesa = n_distinct(ID),.groups = "drop") 
  
  proves_t <- riquesa_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Riquesa[Tractament == "S"]),
      mean_NS = mean(Riquesa[Tractament == "NS"]),
      t_test_result = ifelse(length(unique(Tractament)) == 2, 
                             t.test(Riquesa ~ Tractament)$p.value,
                             NA)
      ,.groups = "drop")
  
  proves_t
## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021     4.33    3.1         0.303 
## 2 Campus-Jun-2021        7.3     9.5         0.190 
## 3 Campus-juliol-2020     3.89    5.92        0.0281
## 4 Campus-maig-2020     NaN       6.55       NA     
## 5 Franqueses-Maig-2021   4       4.6         0.743 
## 6 Sabadell-Juny-2021     3.43    4.8         0.201
ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Riquesa per Zona i Tractaments",
       x = "Mostreig",
       y = "Riquesa")

Generalment es compleix un patró de major riquesa al tractament NS (excepte a Campus-Juliol-2021). Al mostreig del Campus-Juliol_2020 aquesta diferència és significativa (p-valor= 0,0281)

Riquesa d’abelles

# Filtrar les dades d'abelles
  dades_abelles <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]
head(dades_abelles)  
## # A tibble: 6 × 26
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6       159 Campus Juliol         7  2021     10 NS         <NA>         
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
riquesa_x_trampa <- dades_abelles %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Riquesa = n_distinct(ID),.groups = "drop") 
  
  proves_t_abelles <- riquesa_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Riquesa[Tractament == "S"]),
      mean_NS = mean(Riquesa[Tractament == "NS"]),
      t_test_result = ifelse(length(unique(Tractament)) == 2, 
                             t.test(Riquesa ~ Tractament)$p.value,
                             NA)
      ,.groups = "drop")
  
  proves_t_abelles
## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021     3.89    2.1        0.0847 
## 2 Campus-Jun-2021        4.3     5.4        0.338  
## 3 Campus-juliol-2020     2.78    4.46       0.00961
## 4 Campus-maig-2020     NaN       3.84      NA      
## 5 Franqueses-Maig-2021   2       3          0.341  
## 6 Sabadell-Juny-2021     3.14    3.8        0.568
 ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Riquesa d'abelles per Zona i Tractaments",
       x = "Mostreig",
       y = "Riquesa")

Generalment es compleix un patró de major riquesa d’abelles al tractament NS (excepte a Campus-Juliol-2021). Al mostreig del Campus-Juliol_2020 aquesta diferència és significativa (p-valor= 0.0096)

Riquesa de coleòpters

# Filtrar les dades per coleòpters
  dades_coleoptera <- dades_finals[dades_finals$Ordre == "Coleoptera", ]
head(dades_coleoptera)  
## # A tibble: 6 × 26
##   Especimen Zona       Mes    `Mes num`   Any Trampa Tractament Identificacio   
##       <dbl> <chr>      <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>           
## 1        29 Campus     Juliol         7  2021      4 NS         <NA>            
## 2       158 Campus     Juliol         7  2021     10 NS         <NA>            
## 3       183 Campus     Juliol         7  2021      5 S          <NA>            
## 4       101 Franqueses Maig           5  2021      2 NS         <NA>            
## 5       150 Franqueses Maig           5  2021      5 NS         Chiqui bupresti…
## 6       151 Franqueses Maig           5  2021      5 NS         Chiqui bupresti…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
riquesa_x_trampa <- dades_coleoptera %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Riquesa = n_distinct(ID),.groups = "drop") 

  proves_t_coleoptera <- riquesa_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Riquesa[Tractament == "S"]),
      mean_NS = mean(Riquesa[Tractament == "NS"]),
      t_test_result = tryCatch(
      t.test(Riquesa ~ Tractament)$p.value,
      error = function(e) NA
    )
      ,.groups = "drop")
  
  proves_t_coleoptera
## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021     1       1           NA    
## 2 Campus-Jun-2021        2.88    3.7          0.357
## 3 Campus-juliol-2020     1.33    1.25         0.851
## 4 Campus-maig-2020     NaN       2.44        NA    
## 5 Franqueses-Maig-2021   1.5     1            0.500
## 6 Sabadell-Juny-2021     1       1.33         0.423
  ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Riquesa x Zona i Tractaments",
       x = "Mostreig",
       y = "Riquesa")

ERROR: NOT ENOUGH OBSERVATIONS. ALGUNS MOSTREJOS NO TENEN COLEÒPTERS.

Riquesa de sírfids

# Filtrar les dades per sírfids
dades_sirfids <-dades_finals %>% filter(Familia=="Syrphidae")
head(dades_sirfids)  
## # A tibble: 6 × 26
##   Especimen Zona       Mes    `Mes num`   Any Trampa Tractament Identificacio   
##       <dbl> <chr>      <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>           
## 1       226 Campus     Juliol         7  2021      7 S          <NA>            
## 2       265 Campus     Juliol         7  2021     14 NS         <NA>            
## 3       277 Franqueses Maig           5  2021      4 NS         <NA>            
## 4       278 Franqueses Maig           5  2021      4 NS         <NA>            
## 5       279 Franqueses Maig           5  2021      4 NS         <NA>            
## 6       280 Franqueses Maig           5  2021      4 NS         Episyrphus balt…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
riquesa_x_trampa <- dades_sirfids %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Riquesa = n_distinct(ID),.groups = "drop") 
  
  proves_t_sirfids <- riquesa_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Riquesa[Tractament == "S"]),
      mean_NS = mean(Riquesa[Tractament == "NS"]),
      t_test_result = tryCatch(
      t.test(Riquesa ~ Tractament)$p.value,
      error = function(e) NA
    )
      ,.groups = "drop")
  
  proves_t_sirfids
## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl> <lgl>        
## 1 Campus-Juliol-2021        1    1    NA           
## 2 Campus-Jun-2021           1    1    NA           
## 3 Campus-juliol-2020        3    1    NA           
## 4 Campus-maig-2020        NaN    1.18 NA           
## 5 Franqueses-Maig-2021      1    3    NA           
## 6 Sabadell-Juny-2021      NaN    1    NA
ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Riquesa x Zona i Tractaments",
       x = "Mostreig",
       y = "Riquesa")

Hi ha poques dades per obtenir resultats del t-test.

Comparació de l’abundància per mostrejos i tractaments

Abundància total

abundancia_x_trampa <- dades_finals %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Abundancia = n(),.groups = "drop") 
  
  proves_t_abtotal <- abundancia_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Abundancia[Tractament == "S"]),
      mean_NS = mean(Abundancia[Tractament == "NS"]),
      t_test_result = ifelse(length(unique(Tractament)) == 2, 
                             t.test(Abundancia ~ Tractament)$p.value,
                             NA)
      ,.groups = "drop")
  
  proves_t_abtotal
## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021    14.1     4.5          0.130
## 2 Campus-Jun-2021       11.6    18.5          0.113
## 3 Campus-juliol-2020     8.89    8.54         0.904
## 4 Campus-maig-2020     NaN      15.4         NA    
## 5 Franqueses-Maig-2021   4.67    9            0.178
## 6 Sabadell-Juny-2021     7.29    5.8          0.606
 ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Abundancia x Zona i Tractaments",
       x = "Mostreig",
       y = "Abundancia")

No hi ha cap resultat significatiu (p-valor<0.05). No hi ha un patró clar de com varia l’abundància en funció del tractament.

Abundància d’abelles

# Filtrar les dades d'abelles
  dades_abelles <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]
head(dades_abelles)  
## # A tibble: 6 × 26
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6       159 Campus Juliol         7  2021     10 NS         <NA>         
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
abundancia_x_trampa <- dades_abelles %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Abundancia = n(),.groups = "drop") 
  
  proves_t_ababelles <- abundancia_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Abundancia[Tractament == "S"]),
      mean_NS = mean(Abundancia[Tractament == "NS"]),
      t_test_result = ifelse(length(unique(Tractament)) == 2, 
                             t.test(Abundancia ~ Tractament)$p.value,
                             NA)
      ,.groups = "drop")
  
  proves_t_ababelles
## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021    13.7     3.5         0.108 
## 2 Campus-Jun-2021        7.6     8.9         0.589 
## 3 Campus-juliol-2020     7.44    7.08        0.896 
## 4 Campus-maig-2020     NaN       8.68       NA     
## 5 Franqueses-Maig-2021   2       4.4         0.0902
## 6 Sabadell-Juny-2021     7       4.8         0.454
 ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Abundancia x Zona i Tractaments",
       x = "Mostreig",
       y = "Abundancia")

A Franqueses l’abundància d’abelles és significativament superior al tractament NS (p-valor=0.09015). En la resta de mostrejos el patró no és clar.

Abundància de coleòpters

# Filtrar les dades per coleòpters
  dades_coleoptera <- dades_finals[dades_finals$Ordre == "Coleoptera", ]
head(dades_coleoptera)  
## # A tibble: 6 × 26
##   Especimen Zona       Mes    `Mes num`   Any Trampa Tractament Identificacio   
##       <dbl> <chr>      <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>           
## 1        29 Campus     Juliol         7  2021      4 NS         <NA>            
## 2       158 Campus     Juliol         7  2021     10 NS         <NA>            
## 3       183 Campus     Juliol         7  2021      5 S          <NA>            
## 4       101 Franqueses Maig           5  2021      2 NS         <NA>            
## 5       150 Franqueses Maig           5  2021      5 NS         Chiqui bupresti…
## 6       151 Franqueses Maig           5  2021      5 NS         Chiqui bupresti…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
abundancia_x_trampa <- dades_coleoptera %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Abundancia = n(),.groups = "drop") 
  
  proves_t_abcoleoptera <- abundancia_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Abundancia[Tractament == "S"]),
      mean_NS = mean(Abundancia[Tractament == "NS"]),
      t_test_result = tryCatch(
      t.test(Abundancia ~ Tractament)$p.value,
      error = function(e) NA
    )
      ,.groups = "drop")
  
  proves_t_abcoleoptera
## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021     1       1           NA    
## 2 Campus-Jun-2021        4.12    9.2          0.175
## 3 Campus-juliol-2020     1.33    1.25         0.851
## 4 Campus-maig-2020     NaN       6.61        NA    
## 5 Franqueses-Maig-2021   2       4.67         0.496
## 6 Sabadell-Juny-2021     1       1.33         0.423
  ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Abundancia x Zona i Tractaments",
       x = "Mostreig",
       y = "Abundancia")

ERROR: NOT ENOUGH OBSERVATIONS. ALGUNS MOSTREJOS NO TENEN COLEÒPTERS.

Abundància de sírfids

# Filtrar les dades per sírfids
dades_sirfids <-dades_finals %>% filter(Familia=="Syrphidae")
head(dades_sirfids)  
## # A tibble: 6 × 26
##   Especimen Zona       Mes    `Mes num`   Any Trampa Tractament Identificacio   
##       <dbl> <chr>      <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>           
## 1       226 Campus     Juliol         7  2021      7 S          <NA>            
## 2       265 Campus     Juliol         7  2021     14 NS         <NA>            
## 3       277 Franqueses Maig           5  2021      4 NS         <NA>            
## 4       278 Franqueses Maig           5  2021      4 NS         <NA>            
## 5       279 Franqueses Maig           5  2021      4 NS         <NA>            
## 6       280 Franqueses Maig           5  2021      4 NS         Episyrphus balt…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
abundancia_x_trampa <- dades_sirfids %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Abundancia = n_distinct(ID),.groups = "drop") 
  
  proves_t_absirfids <- abundancia_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Abundancia[Tractament == "S"]),
      mean_NS = mean(Abundancia[Tractament == "NS"]),
      t_test_result = tryCatch(
      t.test(Abundancia ~ Tractament)$p.value,
      error = function(e) NA
    )
      ,.groups = "drop")
  
  proves_t_absirfids
## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl> <lgl>        
## 1 Campus-Juliol-2021        1    1    NA           
## 2 Campus-Jun-2021           1    1    NA           
## 3 Campus-juliol-2020        3    1    NA           
## 4 Campus-maig-2020        NaN    1.18 NA           
## 5 Franqueses-Maig-2021      1    3    NA           
## 6 Sabadell-Juny-2021      NaN    1    NA
  ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Abundancia x Zona i Tractaments",
       x = "Mostreig",
       y = "Abundancia")

No hi ha prous dades per obtenir bons resultats de la prova t-student.

DIAGRAMES D’EULER:

Diagrames d’Euler totals per cada zona

# Preparació de les dades: càlcul de l'abundància de cada espècie (ID) a cada mostreig i tractament (NS, S i intersecció)
taula_aux_euler <- as.data.frame.matrix(table(interaction(dades_finals$Zona,dades_finals$ID),dades_finals$Tractament)) %>% mutate(`NS&S`=pmin(NS,S)) 

head(taula_aux_euler)
##                                  NS S NS&S
## Campus.Acmaeodera cylindrica      5 1    1
## Franqueses.Acmaeodera cylindrica  0 0    0
## Sabadell.Acmaeodera cylindrica    0 0    0
## Campus.Amegilla sp.               1 0    0
## Franqueses.Amegilla sp.           0 0    0
## Sabadell.Amegilla sp.             0 0    0
# Reorganitzar la taula auxiliar perquè mostri el nom de les espècies (ID) trobades en cada mostreig i tractament
taula_euler_especies <- taula_aux_euler %>% mutate(Especie = sub("^[^.]+\\.", "", rownames(.)), Zona = sub("\\..*","",rownames(taula_aux_euler)))
taula_euler_especies <- taula_euler_especies %>% group_by(Zona) %>% summarise(  S = list(unique(Especie[ S > 0 ])),
    NS = list(unique(Especie[ NS > 0 ])),
   `NS&S` = list(unique(Especie[ `NS&S` > 0 ]))  )

head(taula_euler_especies)
## # A tibble: 3 × 4
##   Zona       S          NS         `NS&S`    
##   <chr>      <list>     <list>     <list>    
## 1 Campus     <chr [59]> <chr [98]> <chr [39]>
## 2 Franqueses <chr [11]> <chr [20]> <chr [5]> 
## 3 Sabadell   <chr [12]> <chr [16]> <chr [5]>
# Descarregar les llistes d'espècies en un excel
write.xlsx(taula_euler_especies, file.path(directorio, "2.0_Euler_Especies.xlsx"))
# Agrupació per mostrejos (Zona, Mes, Any)
euler_data <- taula_euler_especies %>%
  mutate(NS = map_dbl(NS, length),
         S = map_dbl(S, length),
         `NS&S` = map_dbl(`NS&S`, length))

head(euler_data)
## # A tibble: 3 × 4
##   Zona           S    NS `NS&S`
##   <chr>      <dbl> <dbl>  <dbl>
## 1 Campus        59    98     39
## 2 Franqueses    11    20      5
## 3 Sabadell      12    16      5
# Creació del diagrama d'Euler
apply_euler <- function(row){
  euler(c("NS"=row$NS,"S"=row$S,"NS&S"=row$`NS&S`))
}

euler_data <- euler_data %>% rowwise() %>% mutate(euler_result= list(apply_euler(cur_data())))
# Plot del diagrama d'Euler
generate_euler_plots <- function(euler_result,zona){
  plot(euler_result,main=paste("Zona: ",zona),quantities = list(type = c("percent", "counts")))}

pmap(list(euler_data$euler_result,euler_data$Zona), generate_euler_plots)
## [[1]]

## 
## [[2]]

## 
## [[3]]

Diagrames d’Euler d’abelles (amb dades_abelles) per cada zona

dades_abelles <- dades_abelles[!(dades_abelles$Mes == "maig" & dades_abelles$Any == 2020 & dades_abelles$Zona == "Campus"), ]
# Preparació de les dades: càlcul de l'abundància de cada espècie (ID) a cada mostreig i tractament (NS, S i intersecció)
taula_aux_eulerabelles <- as.data.frame.matrix(table(interaction(dades_abelles$Zona,dades_abelles$ID),dades_abelles$Tractament)) %>% mutate(`NS&S`=pmin(NS,S)) 

head(taula_aux_eulerabelles)
##                            NS S NS&S
## Campus.Amegilla sp.         1 0    0
## Franqueses.Amegilla sp.     0 0    0
## Sabadell.Amegilla sp.       0 0    0
## Campus.Andrena cinerea      0 0    0
## Franqueses.Andrena cinerea  1 0    0
## Sabadell.Andrena cinerea    0 0    0
# Reorganitzar la taula auxiliar perquè mostri el nom de les espècies (ID) trobades en cada mostreig i tractament
taula_euler_abelles <- taula_aux_eulerabelles %>% mutate(Especie = sub("^[^.]+\\.", "", rownames(.)), Zona = sub("\\..*","",rownames(taula_aux_eulerabelles)))
taula_euler_abelles <- taula_euler_abelles %>% group_by(Zona) %>% summarise(  S = list(unique(Especie[ S > 0 ])),
    NS = list(unique(Especie[ NS > 0 ])),
   `NS&S` = list(unique(Especie[ `NS&S` > 0 ]))  )

head(taula_euler_abelles)
## # A tibble: 3 × 4
##   Zona       S          NS         `NS&S`    
##   <chr>      <list>     <list>     <list>    
## 1 Campus     <chr [35]> <chr [47]> <chr [21]>
## 2 Franqueses <chr [5]>  <chr [12]> <chr [3]> 
## 3 Sabadell   <chr [10]> <chr [12]> <chr [5]>
# Descarregar les llistes d'espècies en un excel
write.xlsx(taula_euler_abelles, file.path(directorio, "Euler_Abelles.xlsx"))
# Agrupació per mostrejos (Zona, Mes, Any)
euler_data_abelles <- taula_euler_abelles %>%
  mutate(NS = map_dbl(NS, length),
         S = map_dbl(S, length),
         `NS&S` = map_dbl(`NS&S`, length))

head(euler_data_abelles)
## # A tibble: 3 × 4
##   Zona           S    NS `NS&S`
##   <chr>      <dbl> <dbl>  <dbl>
## 1 Campus        35    47     21
## 2 Franqueses     5    12      3
## 3 Sabadell      10    12      5
# Creació del diagrama d'Euler
apply_euler <- function(row){
  euler(c("NS"=row$NS,"S"=row$S,"NS&S"=row$`NS&S`))
}

euler_data_abelles <- euler_data_abelles %>% rowwise() %>% mutate(euler_result= list(apply_euler(cur_data())))
# Plot del diagrama d'Euler
generate_euler_plots <- function(euler_result,zona){
  plot(euler_result,main=paste("Zona: ",zona),quantities = list(type = c("percent", "counts")))}

pmap(list(euler_data_abelles$euler_result,euler_data_abelles$Zona), generate_euler_plots)
## [[1]]

## 
## [[2]]

## 
## [[3]]

ANÀLISIS DE COMPONENTS PRINCIPALS (PCA)

PCA total

# Preparar dades: taula de dades Shannon per mostreig i tractament (MT). Normalitzar les dades: nombre d'individus entre el nombre de pantraps del mostreig 
pca_data <- data_shannon_MT

# Calcular el nombre de trampes per cada mostreig i tractament
n_trampes <- dades_finals %>% group_by(Mostreig, Tractament) %>% summarise(numtrampes=n_distinct(Trampa)) %>% ungroup() %>%  mutate(row=paste(Mostreig,Tractament,sep="_")) %>% dplyr::select(-Mostreig,-Tractament)
## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.
# Unir les dues taules relacionant les files del mateix mostreig i tractament
pca_data$row <- rownames(pca_data)
pca_data <- pca_data  %>% inner_join(n_trampes, by=join_by(row))

# Dividir les abundàncies de cada ID entre el nombre de trampes de cada mostreig i tractament
num_columns <-ncol(pca_data) 
pca_data[, 1:(num_columns - 2)] <- pca_data[, 1:(num_columns - 2)] / pca_data[, num_columns]

# Eliminar l'última columna de la taula (numtrampes), ja no cal per fer la PCA
pca_data <- pca_data %>% dplyr::select(-numtrampes)

pca_data<- pca_data %>%
  mutate(Tractament = str_split(row, "_", simplify = TRUE)[, 2])
# Fer la PCA
pca_comps <- prcomp(pca_data %>% dplyr::select(-row,-Tractament))
# Mostrar la PCA
ggbiplot(pca_comps, labels = pca_data$row, groups = pca_data$Tractament, 
         ellipse = TRUE, var.axes = FALSE) +
  theme_minimal() +
  theme(text = element_text(size = 20))

PCA abelles

# Preparar dades: taula de dades Shannon d'abelles. Normalitzar les dades: nombre d'individus entre el nombre de pantraps del mostreig 
pca_data_abelles <- data_shannon_abelles

# Calcular el nombre de trampes per cada mostreig i tractament
n_trampes_abelles <- dades_finals %>%
  filter(Ordre == 'HymenopteraAbella') %>%
  group_by(Mostreig, Tractament) %>%
  summarise(numtrampes = n_distinct(Trampa)) %>%
  ungroup() %>%
  mutate(row = paste(Mostreig, Tractament, sep = "_")) %>%
  dplyr::select(-Mostreig, -Tractament)
## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.
# Unir les dues taules relacionant les files del mateix mostreig i tractament
pca_data_abelles$row <- rownames(pca_data_abelles)
pca_data_abelles <- pca_data_abelles  %>% inner_join(n_trampes_abelles, by=join_by(row))

# Dividir les abundàncies de cada ID entre el nombre de trampes de cada mostreig i tractament
num_columnsabe <-ncol(pca_data_abelles) 
pca_data_abelles[, 1:(num_columnsabe - 2)] <- pca_data_abelles[, 1:(num_columnsabe - 2)] / pca_data_abelles[, num_columnsabe]

# Eliminar l'última columna de la taula (numtrampes), ja no cal per fer la PCA
pca_data_abelles <- pca_data_abelles %>% dplyr::select(-numtrampes)

pca_data_abelles<- pca_data_abelles %>%
  mutate(Tractament = str_split(row, "_", simplify = TRUE)[, 2])

View(pca_data_abelles)
# Fer la PCA
pca_abelles <- prcomp(pca_data_abelles %>% dplyr::select(-row,-Tractament))
summary(pca_abelles)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     1.9981 1.1729 0.61490 0.52897 0.48323 0.38196 0.31303
## Proportion of Variance 0.6025 0.2076 0.05706 0.04223 0.03524 0.02202 0.01479
## Cumulative Proportion  0.6025 0.8101 0.86717 0.90940 0.94464 0.96665 0.98144
##                           PC8     PC9    PC10      PC11
## Standard deviation     0.2482 0.19318 0.15503 2.126e-16
## Proportion of Variance 0.0093 0.00563 0.00363 0.000e+00
## Cumulative Proportion  0.9907 0.99637 1.00000 1.000e+00
# Mostrar la PCA
ggbiplot(pca_abelles, labels = pca_data_abelles$row, groups = pca_data_abelles$Tractament, 
         ellipse = TRUE, var.axes = FALSE) +
  theme_minimal() +
  theme(text = element_text(size = 20))

PCA abelles Campus (per trampes. No sé ben bé què he fet)

# Preparar dades: taula de dades Shannon d'abelles per mostreig, tractament i trampa (dades NO normalitzades)
pca_data_campus <- data_shannon_MTTR %>% 
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  filter(Zona == 'Campus') %>% 
  mutate(Tractament = str_extract(rownames(.), "(?<=_)[^_]+(?=_)"))

#View(pca_data_campus)

# Fer la PCA
pca_campus <- prcomp(pca_data_campus %>% dplyr::select(-Zona,-Tractament))

# Mostrar la PCA
ggbiplot(pca_campus, labels = pca_data_campus$row, groups = pca_data_campus$Tractament, 
         ellipse = TRUE, var.axes = FALSE) +
  theme_minimal() +
  theme(text = element_text(size = 20))

NMDS, ANÀLISIS DE SIMILITUD (ANOSIM), ADONIS, I ESPÈCIES INDICADORES

NMDS total

# Preparar dades
NMDS_data <- data_shannon_MT

# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata <- NMDS_data %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Tractament_metadata <- Tractament_metadata %>% mutate(color=ifelse(Tractament=="NS","green","orange"))

Tractament_metadata
##                         Tractament  color
## Campus-juliol-2020_NS           NS  green
## Campus-juliol-2020_S             S orange
## Campus-Juliol-2021_NS           NS  green
## Campus-Juliol-2021_S             S orange
## Campus-Jun-2021_NS              NS  green
## Campus-Jun-2021_S                S orange
## Campus-maig-2020_NS             NS  green
## Franqueses-Maig-2021_NS         NS  green
## Franqueses-Maig-2021_S           S orange
## Sabadell-Juny-2021_NS           NS  green
## Sabadell-Juny-2021_S             S orange
# NMDS Total
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)

set.seed(7)

# k=2: nombre de dimensions
NMDS <- metaMDS(NMDS_data, k=2)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1574807 
## Run 1 stress 0.1574807 
## ... Procrustes: rmse 1.586753e-06  max resid 1.972763e-06 
## ... Similar to previous best
## Run 2 stress 0.1332528 
## ... New best solution
## ... Procrustes: rmse 0.2375972  max resid 0.4227959 
## Run 3 stress 0.1332526 
## ... New best solution
## ... Procrustes: rmse 0.0002214546  max resid 0.0003933833 
## ... Similar to previous best
## Run 4 stress 0.1271718 
## ... New best solution
## ... Procrustes: rmse 0.07482385  max resid 0.1553409 
## Run 5 stress 0.1927377 
## Run 6 stress 0.1332526 
## Run 7 stress 0.1576285 
## Run 8 stress 0.1271718 
## ... Procrustes: rmse 6.804317e-05  max resid 0.0001405358 
## ... Similar to previous best
## Run 9 stress 0.3407341 
## Run 10 stress 0.1271718 
## ... Procrustes: rmse 1.536335e-05  max resid 3.205959e-05 
## ... Similar to previous best
## Run 11 stress 0.1332526 
## Run 12 stress 0.1576277 
## Run 13 stress 0.1271718 
## ... Procrustes: rmse 4.444601e-06  max resid 9.072617e-06 
## ... Similar to previous best
## Run 14 stress 0.1647707 
## Run 15 stress 0.1271718 
## ... Procrustes: rmse 1.498339e-05  max resid 3.060948e-05 
## ... Similar to previous best
## Run 16 stress 0.1574807 
## Run 17 stress 0.1332526 
## Run 18 stress 0.1576295 
## Run 19 stress 0.1332525 
## Run 20 stress 0.1927396 
## *** Best solution repeated 4 times
# Plot NMDS total
ordiplot(NMDS,type="n")
ordihull(NMDS,groups=Tractament_metadata$Tractament,draw="polygon",col="grey90",label=F)
orditorp(NMDS,display="species",col="black",air=0.01)
orditorp(NMDS,display="sites",col=Tractament_metadata$color,
   air=0.01,cex=1.25)

# Save the last plot
g <- recordPlot()

# Replay and save as PNG
width <- 17 # Width in inches
height <- 14  # Height in inches
res <- 500   # Resolution in dpi

# Replay and save as PNG with higher resolution and larger size
png("Tot.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()
## quartz_off_screen 
##                 2
Descripción del gráfico
Descripción del gráfico

Anàlisi de similitud (ANOSIM) Total

ano = anosim(NMDS_data, Tractament_metadata$Tractament, distance = "bray", permutations = 9999)
ano
## 
## Call:
## anosim(x = NMDS_data, grouping = Tractament_metadata$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: -0.06667 
##       Significance: 0.721 
## 
## Permutation: free
## Number of permutations: 9999

No signiticatiu (p valor > 0.05): la composició de les comunitats no és significativament diferent entre els tractaments S i NS. També s’observa al plot de NMDS: els polígons que corresponen als dos tractaments estan superposats. Això probablement sigui per haver fet la prova amb totes les zones de mostreig alhora.

ADONIS total

# Extreure de la taula el tractament de cada mostreig: una nova columna
Zona_Tractament_metadata_TOTAL <- NMDS_data %>%
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Zona, Tractament)

Zona_Tractament_metadata_TOTAL
##                               Zona Tractament
## Campus-juliol-2020_NS       Campus         NS
## Campus-juliol-2020_S        Campus          S
## Campus-Juliol-2021_NS       Campus         NS
## Campus-Juliol-2021_S        Campus          S
## Campus-Jun-2021_NS          Campus         NS
## Campus-Jun-2021_S           Campus          S
## Campus-maig-2020_NS         Campus         NS
## Franqueses-Maig-2021_NS Franqueses         NS
## Franqueses-Maig-2021_S  Franqueses          S
## Sabadell-Juny-2021_NS     Sabadell         NS
## Sabadell-Juny-2021_S      Sabadell          S
NMDStot.dist <- vegdist(wisconsin(sqrt(NMDS_data)), k=2)

adonis(formula = NMDStot.dist~Tractament/Zona, data = Zona_Tractament_metadata_TOTAL, permutations = 1000)$aov.tab
## 'adonis' will be deprecated: use 'adonis2' instead
## Permutation: free
## Number of permutations: 1000
## 
## Terms added sequentially (first to last)
## 
##                 Df SumsOfSqs MeanSqs F.Model      R2  Pr(>F)  
## Tractament       1    0.3611 0.36111  1.1373 0.10191 0.25874  
## Tractament:Zona  4    1.5947 0.39868  1.2556 0.45004 0.03996 *
## Residuals        5    1.5877 0.31753         0.44805          
## Total           10    3.5435                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No hi ha diferències significatives entre tractaments (p-valor = 0.25874), però la interacció entre zona i tractament explica un 45% de la variabilitat, i és un resultat significatiu (p-valor=0.03996).

Espècies indicadores dels tractaments S i NS total

indicadores = multipatt(NMDS_data, Tractament_metadata$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 137
##  Selected number of species: 0 
##  Number of species associated to 1 group: 0 
## 
##  List of species associated to each combination: 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No hi ha espècies indicadores significatives de cap dels dos tractaments.

summary(indicadores, alpha = 0.2)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.2
## 
##  Total number of species: 137
##  Selected number of species: 6 
##  Number of species associated to 1 group: 6 
## 
##  List of species associated to each combination: 
## 
##  Group NS  #sps.  3 
##                    stat p.value
## Anthaxia sp.      0.422   0.182
## Eucera elongatula 0.422   0.181
## Ceratina parvula  0.391   0.178
## 
##  Group S  #sps.  3 
##                             stat p.value
## Vespula germanica          0.500   0.186
## Lasioglossum glabriusculum 0.430   0.184
## Lasioglossum griseolum     0.385   0.184
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Si augmentem el nivell de significació (alfa) a 0,2, hi ha 3 espècies indicadores de NS (Anthaxia sp., Eucera elongatula, Ceratina parvula), i 3 espècies indicadores de S (Vespula germanica, Lasioglossum gabriusculum i Lasioglossum griseolum). Però realment no son significatives.

NMDS Campus

# Filtrar les dades del Campus i preparar les dades
NMDS_data_campus <- NMDS_data %>%
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  filter(Zona=="Campus") %>% dplyr::select(-Zona)

# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata <- NMDS_data_campus %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Tractament_metadata <- Tractament_metadata %>% mutate(color=ifelse(Tractament=="NS","green","orange"))

# NMDS_data_campus
# NMDS Campus
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)

set.seed(7)

# k=2: nombre de dimensions
NMDScampus <- metaMDS(NMDS_data_campus, k=2)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.03999251 
## Run 1 stress 0.03999252 
## ... Procrustes: rmse 3.432746e-05  max resid 4.816299e-05 
## ... Similar to previous best
## Run 2 stress 0.03999251 
## ... New best solution
## ... Procrustes: rmse 1.138143e-05  max resid 1.575934e-05 
## ... Similar to previous best
## Run 3 stress 0.03999257 
## ... Procrustes: rmse 0.0001319502  max resid 0.0001794008 
## ... Similar to previous best
## Run 4 stress 0.2570537 
## Run 5 stress 0.03999251 
## ... Procrustes: rmse 1.730368e-05  max resid 2.266611e-05 
## ... Similar to previous best
## Run 6 stress 0.04217697 
## Run 7 stress 0.042177 
## Run 8 stress 0.2570032 
## Run 9 stress 0.04581483 
## Run 10 stress 0.04581482 
## Run 11 stress 0.2570537 
## Run 12 stress 0.04318443 
## Run 13 stress 0.2570537 
## Run 14 stress 0.04318432 
## Run 15 stress 0.0432258 
## Run 16 stress 0.03999251 
## ... Procrustes: rmse 1.482389e-05  max resid 2.245167e-05 
## ... Similar to previous best
## Run 17 stress 0.293306 
## Run 18 stress 0.04581486 
## Run 19 stress 0.04625739 
## Run 20 stress 0.1960094 
## *** Best solution repeated 4 times
# Plot NMDS Campus
ordiplot(NMDScampus,type="n")
ordihull(NMDScampus,groups=Tractament_metadata$Tractament,draw="polygon",col="grey90",label=F)
orditorp(NMDScampus,display="species",col="black",air=0.01)
orditorp(NMDScampus,display="sites",col=Tractament_metadata$color,
   air=0.01,cex=1.25)

# Save the last plot
g <- recordPlot()

# Replay and save as PNG
width <- 17 # Width in inches
height <- 14  # Height in inches
res <- 500   # Resolution in dpi

# Replay and save as PNG with higher resolution and larger size
png("output_plot.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()
Descripción del gráfico
Descripción del gráfico

Anàlisi de similitud (ANOSIM) Campus

ano = anosim(NMDS_data_campus, Tractament_metadata$Tractament, distance = "bray", permutations = 9999)
## 'nperm' >= set of all permutations: complete enumeration.
## Set of permutations < 'minperm'. Generating entire set.
ano
## 
## Call:
## anosim(x = NMDS_data_campus, grouping = Tractament_metadata$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: -0.03704 
##       Significance: 0.42857 
## 
## Permutation: free
## Number of permutations: 5039

No signiticatiu (p valor = 0.42857): la composició de les comunitats de S i NS al Campus no és significativament diferent entre els dos tractaments (S i NS).

ADONIS Campus

# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata_CAMPUS <- NMDS_data_campus %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

Tractament_metadata_CAMPUS
##                       Tractament
## Campus-juliol-2020_NS         NS
## Campus-juliol-2020_S           S
## Campus-Juliol-2021_NS         NS
## Campus-Juliol-2021_S           S
## Campus-Jun-2021_NS            NS
## Campus-Jun-2021_S              S
## Campus-maig-2020_NS           NS
NMDScamp.dist <- vegdist(wisconsin(sqrt(NMDS_data_campus)), k=2)

adonis(formula = NMDScamp.dist~Tractament, data = Tractament_metadata_CAMPUS, permutations = 1000)$aov.tab
## 'adonis' will be deprecated: use 'adonis2' instead
## Set of permutations < 'minperm'. Generating entire set.
## Permutation: free
## Number of permutations: 5039
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## Tractament  1   0.27747 0.27747 0.87044 0.14828 0.5924
## Residuals   5   1.59387 0.31877         0.85172       
## Total       6   1.87134                 1.00000

El tractament només explica un 14,8% de la variabilitat, i el resultat no és significatiu.

Espècies indicadores dels tractaments S i NS Campus

indicadores = multipatt(NMDS_data_campus, Tractament_metadata$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 137
##  Selected number of species: 0 
##  Number of species associated to 1 group: 0 
## 
##  List of species associated to each combination: 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No hi ha espècies indicadores significatives de cap dels dos tractaments.

summary(indicadores, alpha = 0.2)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.2
## 
##  Total number of species: 137
##  Selected number of species: 5 
##  Number of species associated to 1 group: 5 
## 
##  List of species associated to each combination: 
## 
##  Group NS  #sps.  3 
##                     stat p.value
## Halictus scabiosae 0.775   0.149
## Anthaxia sp.       0.542   0.147
## Ceratina parvula   0.498   0.147
## 
##  Group S  #sps.  2 
##                             stat p.value
## Lasioglossum glabriusculum 0.731   0.118
## Vespula germanica          0.707   0.145
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Si augmentem el nivell de significació (alfa) a 0,2, hi ha 3 espècies indicadores de NS (Halictus scabiosae, Anthaxia sp., i Ceratina parvula), i 2 espècies indicadores de S (Lasioglossum glabriusculum i Vespula germanica). Però no son resultats significatius.

NMDS d’abelles totals (amb data_shannon_abelles)

# Crear metadades del tractament
NMDSabelles_data <- data_shannon_abelles

# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata <- NMDSabelles_data %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Tractament_metadata <- Tractament_metadata %>% mutate(color=ifelse(Tractament=="NS","green","orange"))

Tractament_metadata
##                         Tractament  color
## Campus-juliol-2020_NS           NS  green
## Campus-juliol-2020_S             S orange
## Campus-Juliol-2021_NS           NS  green
## Campus-Juliol-2021_S             S orange
## Campus-Jun-2021_NS              NS  green
## Campus-Jun-2021_S                S orange
## Campus-maig-2020_NS             NS  green
## Franqueses-Maig-2021_NS         NS  green
## Franqueses-Maig-2021_S           S orange
## Sabadell-Juny-2021_NS           NS  green
## Sabadell-Juny-2021_S             S orange
# NMDS abelles
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)

set.seed(7)

# k=2: nombre de dimensions
NMDSabe <- metaMDS(NMDSabelles_data, k=2)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1466713 
## Run 1 stress 0.1762197 
## Run 2 stress 0.1764663 
## Run 3 stress 0.1430526 
## ... New best solution
## ... Procrustes: rmse 0.1296798  max resid 0.3307691 
## Run 4 stress 0.1660287 
## Run 5 stress 0.1350214 
## ... New best solution
## ... Procrustes: rmse 0.08009816  max resid 0.2022004 
## Run 6 stress 0.1430526 
## Run 7 stress 0.1660287 
## Run 8 stress 0.1350213 
## ... New best solution
## ... Procrustes: rmse 0.000511821  max resid 0.0009819792 
## ... Similar to previous best
## Run 9 stress 0.3407333 
## Run 10 stress 0.1350213 
## ... Procrustes: rmse 0.0003994713  max resid 0.0007627691 
## ... Similar to previous best
## Run 11 stress 0.1466713 
## Run 12 stress 0.1350216 
## ... Procrustes: rmse 0.0005481505  max resid 0.001090682 
## ... Similar to previous best
## Run 13 stress 0.1430527 
## Run 14 stress 0.1487731 
## Run 15 stress 0.1430526 
## Run 16 stress 0.1660287 
## Run 17 stress 0.1487731 
## Run 18 stress 0.211699 
## Run 19 stress 0.1466713 
## Run 20 stress 0.1350213 
## ... Procrustes: rmse 0.0002187493  max resid 0.000402306 
## ... Similar to previous best
## *** Best solution repeated 4 times
NMDSabe
## 
## Call:
## metaMDS(comm = NMDSabelles_data, k = 2) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(NMDSabelles_data)) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.1350213 
## Stress type 1, weak ties
## Best solution was repeated 4 times in 20 tries
## The best solution was from try 8 (random start)
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(NMDSabelles_data))'
# Plot NMDS abelles
ordiplot(NMDSabe,type="n")
ordihull(NMDSabe,groups=Tractament_metadata$Tractament,draw="polygon",col="grey90",label=F)
orditorp(NMDSabe,display="species",col="black",air=0.01)
orditorp(NMDSabe,display="sites",col=Tractament_metadata$color,
   air=0.01,cex=1.25)

# Save the last plot
g <- recordPlot()

# Replay and save as PNG
width <- 17 # Width in inches
height <- 14  # Height in inches
res <- 500   # Resolution in dpi

# Replay and save as PNG with higher resolution and larger size
png("AbellesTot.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()
Descripción del gráfico
Descripción del gráfico

Anàlisi de similitud (ANOSIM) abelles totals

ano = anosim(NMDSabelles_data, Tractament_metadata$Tractament, distance = "bray", permutations = 9999)
ano
## 
## Call:
## anosim(x = NMDSabelles_data, grouping = Tractament_metadata$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: -0.04267 
##       Significance: 0.6117 
## 
## Permutation: free
## Number of permutations: 9999

No signiticatiu (p valor > 0.6037): la composició de les comunitats és pràcticament idèntica (valor de l’estadístic R proper a 0, i fins i tot negatiu). A més, el resultat no és significatiu. També s’observa al plot de NMDS: els polígons que corresponen als dos tractaments estan superposats. Això probablement sigui per haver fet la prova amb totes les zones de mostreig alhora.

ADONIS abelles totals

# Extreure de la taula el tractament de cada mostreig: una nova columna
Zona_Tractament_metadata <- NMDSabelles_data %>%
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Zona, Tractament)

Zona_Tractament_metadata
##                               Zona Tractament
## Campus-juliol-2020_NS       Campus         NS
## Campus-juliol-2020_S        Campus          S
## Campus-Juliol-2021_NS       Campus         NS
## Campus-Juliol-2021_S        Campus          S
## Campus-Jun-2021_NS          Campus         NS
## Campus-Jun-2021_S           Campus          S
## Campus-maig-2020_NS         Campus         NS
## Franqueses-Maig-2021_NS Franqueses         NS
## Franqueses-Maig-2021_S  Franqueses          S
## Sabadell-Juny-2021_NS     Sabadell         NS
## Sabadell-Juny-2021_S      Sabadell          S
NMDSabe.dist <- vegdist(wisconsin(sqrt(NMDSabelles_data)), k=2)

adonis(formula = NMDSabe.dist~Tractament/Zona, data = Zona_Tractament_metadata, permutations = 1000)$aov.tab
## 'adonis' will be deprecated: use 'adonis2' instead
## Permutation: free
## Number of permutations: 1000
## 
## Terms added sequentially (first to last)
## 
##                 Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## Tractament       1    0.3748 0.37477  1.3134 0.11721 0.1299  
## Tractament:Zona  4    1.3958 0.34896  1.2229 0.43656 0.0979 .
## Residuals        5    1.4267 0.28534         0.44622         
## Total           10    3.1973                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El tractament només explica un 11,72% de la variabilitat observada entre les poblacions de S i NS, i no és un resultat significatiu. La interacció entre el tractament i la zona explica un 43,65% de la variabilitat observada, i és un resultat marginalment significatiu (0,0989, < 0.1).

Espècies d’abelles indicadores dels tractaments S i NS

indicadores = multipatt(NMDSabelles_data, Tractament_metadata$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 80
##  Selected number of species: 0 
##  Number of species associated to 1 group: 0 
## 
##  List of species associated to each combination: 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No hi ha espècies indicadores significatives.

summary(indicadores, alpha = 0.2)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.2
## 
##  Total number of species: 80
##  Selected number of species: 3 
##  Number of species associated to 1 group: 3 
## 
##  List of species associated to each combination: 
## 
##  Group NS  #sps.  1 
##                    stat p.value
## Eucera elongatula 0.422   0.181
## 
##  Group S  #sps.  2 
##                             stat p.value
## Lasioglossum glabriusculum 0.430   0.184
## Lasioglossum griseolum     0.385   0.184
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Si augmentem el nivell de significació (alfa) a 0.2, hi ha 1 espècie indicadora a NS (Eucera elongatula) i 2 espècies indicadores a S (L. glabriusculum i L. griseolum). Però no son resultats significatius.

NMDS abelles al Campus:

# Filtrar les dades del Campus i preparar les dades
NMDSabelles_data_campus <- NMDSabelles_data %>%
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  filter(Zona=="Campus") %>% dplyr::select(-Zona)

# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata <- NMDSabelles_data_campus %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Tractament_metadata <- Tractament_metadata %>% mutate(color=ifelse(Tractament=="NS","green","orange"))
# NMDS abelles Campus
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)

set.seed(7)

# k=2: nombre de dimensions
NMDSabe_Campus <- metaMDS(NMDSabelles_data_campus, k=2)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.05473171 
## Run 1 stress 0.05473171 
## ... Procrustes: rmse 1.536864e-06  max resid 2.917548e-06 
## ... Similar to previous best
## Run 2 stress 0.09163904 
## Run 3 stress 0.05587716 
## Run 4 stress 0.2601494 
## Run 5 stress 0.05473171 
## ... Procrustes: rmse 3.481031e-06  max resid 6.748059e-06 
## ... Similar to previous best
## Run 6 stress 0.05473171 
## ... New best solution
## ... Procrustes: rmse 1.075571e-06  max resid 1.942485e-06 
## ... Similar to previous best
## Run 7 stress 0.2143136 
## Run 8 stress 0.05473171 
## ... Procrustes: rmse 3.503063e-06  max resid 6.969015e-06 
## ... Similar to previous best
## Run 9 stress 0.05587716 
## Run 10 stress 0.1001423 
## Run 11 stress 0.05587716 
## Run 12 stress 0.05473171 
## ... Procrustes: rmse 3.712786e-06  max resid 7.351559e-06 
## ... Similar to previous best
## Run 13 stress 0.05587716 
## Run 14 stress 0.05473171 
## ... Procrustes: rmse 3.16519e-06  max resid 6.1088e-06 
## ... Similar to previous best
## Run 15 stress 0.2294233 
## Run 16 stress 0.1001423 
## Run 17 stress 0.293306 
## Run 18 stress 0.05587716 
## Run 19 stress 0.1001423 
## Run 20 stress 0.1964881 
## *** Best solution repeated 4 times
# Plot NMDS
ordiplot(NMDSabe_Campus,type="n")
ordihull(NMDSabe_Campus,groups=Tractament_metadata$Tractament,draw="polygon",col="grey90",label=F)
orditorp(NMDSabe_Campus,display="species",col="black",air=0.01)
orditorp(NMDSabe_Campus,display="sites",col=Tractament_metadata$color,
   air=0.01,cex=1.25)

# Save the last plot
g <- recordPlot()

# Replay and save as PNG
width <- 17 # Width in inches
height <- 14  # Height in inches
res <- 500   # Resolution in dpi

# Replay and save as PNG with higher resolution and larger size
png("AbellesCampus.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()
Descripción del gráfico
Descripción del gráfico

Anàlisi de similitud (ANOSIM) abelles Campus

ano = anosim(NMDSabelles_data_campus, Tractament_metadata$Tractament, distance = "bray", permutations = 9999)
## 'nperm' >= set of all permutations: complete enumeration.
## Set of permutations < 'minperm'. Generating entire set.
ano
## 
## Call:
## anosim(x = NMDSabelles_data_campus, grouping = Tractament_metadata$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.03704 
##       Significance: 0.34286 
## 
## Permutation: free
## Number of permutations: 5039

La composició de les comunitats entre els dos tractaments és molt similar (valor de l’estadístic d’R molt proper a 0), i el resultat no és significatiu.

ADONIS abelles Campus

# Extreure de la taula el tractament de cada mostreig: una nova columna
ADONIS1_Tractament_metadata <- NMDSabelles_data_campus %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

ADONIS1_Tractament_metadata
##                       Tractament
## Campus-juliol-2020_NS         NS
## Campus-juliol-2020_S           S
## Campus-Juliol-2021_NS         NS
## Campus-Juliol-2021_S           S
## Campus-Jun-2021_NS            NS
## Campus-Jun-2021_S              S
## Campus-maig-2020_NS           NS
NMDSabe.dist <- vegdist(wisconsin(sqrt(NMDSabelles_data_campus)), k=2)

adonis(formula = NMDSabe.dist~Tractament, data = ADONIS1_Tractament_metadata, permutations = 1000)$aov.tab
## 'adonis' will be deprecated: use 'adonis2' instead
## Set of permutations < 'minperm'. Generating entire set.
## Permutation: free
## Number of permutations: 5039
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## Tractament  1   0.26326 0.26326  0.9175 0.15505 0.4905
## Residuals   5   1.43466 0.28693         0.84495       
## Total       6   1.69792                 1.00000

El factor Tractament explica un 15.5% de la variabilitat d’abelles al Campus entre els tractaments S i NS, però aquesta explicació no és significativa.

Espècies d’abelles indicadores dels tractaments S i NS al Campus

indicadores = multipatt(NMDSabelles_data_campus, Tractament_metadata$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores, alpha = 0.2)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.2
## 
##  Total number of species: 80
##  Selected number of species: 2 
##  Number of species associated to 1 group: 2 
## 
##  List of species associated to each combination: 
## 
##  Group NS  #sps.  1 
##                     stat p.value
## Halictus scabiosae 0.775   0.149
## 
##  Group S  #sps.  1 
##                             stat p.value
## Lasioglossum glabriusculum 0.731   0.118
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No hi ha espècies indicadores significatives. Si augmenta el nivell de significació alfa = 0.2, sí que apareixen Halictus scabiosae com a indicador del tractament NO segat, i Lasioglossum glabriusculum com a indicador del tractament segat.

NMDS amb presència-absència (no abundàncies) ABELLES CAMPUS

NMDS_data_abellescampus_presabs<-NMDSabelles_data_campus
NMDS_data_abellescampus_presabs[NMDS_data_abellescampus_presabs>0]=1
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)

set.seed(7)

# k=2: nombre de dimensions
NMDS_presabs <- metaMDS(NMDS_data_abellescampus_presabs, k=2)

ordiplot(NMDS_presabs,type="n")
ordihull(NMDS_presabs,groups=Tractament_metadata$Tractament,draw="polygon",col="grey90",label=F)
orditorp(NMDS_presabs,display="species",col="black",air=0.01)
orditorp(NMDS_presabs,display="sites",col=Tractament_metadata$color,
   air=0.01,cex=1.25)

# Save the last plot
g <- recordPlot()

# Replay and save as PNG
width <- 17 # Width in inches
height <- 14  # Height in inches
res <- 500   # Resolution in dpi

# Replay and save as PNG with higher resolution and larger size
png("Presenciabsencia.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()
Descripción del gráfico
Descripción del gráfico

Anàlisi de similitud ANOSIM presència-absència

ano = anosim(NMDS_data_abellescampus_presabs, Tractament_metadata$Tractament, distance = "bray", permutations = 9999)
## 'nperm' >= set of all permutations: complete enumeration.
## Set of permutations < 'minperm'. Generating entire set.
ano
## 
## Call:
## anosim(x = NMDS_data_abellescampus_presabs, grouping = Tractament_metadata$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: -0.1019 
##       Significance: 0.62857 
## 
## Permutation: free
## Number of permutations: 5039

Les composicions de les comunitats dels dos tractaments son molt similars, i el resultat no és significatiu.

ADONIS presència-absència

# Extreure de la taula el tractament de cada mostreig: una nova columna
ADONIS_Tractament_metadata <- NMDS_data_abellescampus_presabs %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

ADONIS_Tractament_metadata
##                       Tractament
## Campus-juliol-2020_NS         NS
## Campus-juliol-2020_S           S
## Campus-Juliol-2021_NS         NS
## Campus-Juliol-2021_S           S
## Campus-Jun-2021_NS            NS
## Campus-Jun-2021_S              S
## Campus-maig-2020_NS           NS
NMDSabe.dist <- vegdist(wisconsin(sqrt(NMDS_data_abellescampus_presabs)), k=2)

adonis(formula = NMDSabe.dist~Tractament, data = ADONIS_Tractament_metadata, permutations = 1000)$aov.tab
## 'adonis' will be deprecated: use 'adonis2' instead
## Set of permutations < 'minperm'. Generating entire set.
## Permutation: free
## Number of permutations: 5039
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## Tractament  1   0.23205 0.23205 0.97825 0.16363 0.3716
## Residuals   5   1.18603 0.23721         0.83637       
## Total       6   1.41808                 1.00000

El tractament explica un 16% de la variabilitat entre les comunitats, però el resultat no és significatiu.

Espècies indicadores presència-absència

indicadores = multipatt(NMDS_data_abellescampus_presabs, Tractament_metadata$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores, alpha = 0.2)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.2
## 
##  Total number of species: 80
##  Selected number of species: 1 
##  Number of species associated to 1 group: 1 
## 
##  List of species associated to each combination: 
## 
##  Group NS  #sps.  1 
##                     stat p.value
## Halictus scabiosae 0.775   0.137
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Només Halictus scabiosae com a indicadora del tractament NS, però el resultat no és significatiu (nivell de significació alfa = 0.2)

MODELS LINEALS MIXTOS D’ANÀLISIS DESCRIPTIVES:

Abundància de Lasioglossums

# Crear matriu amb només les columnes de Lasioglossums

noms_columnes <- colnames(metadata_shannon)
genere <- "Lasioglossum"

columnes_lasioglossum <- noms_columnes[grep(paste0("^", genere), noms_columnes)]

dades_lasioglossum <- metadata_shannon[, columnes_lasioglossum]

# Agrupar totes les columnes en una sola (abundància total de Lasioglossums en cada trampa)
dades_lasioglossum$Abundancia_Total <- rowSums(dades_lasioglossum)
dades_lasioglossum_ok <- data.frame(Abundancia_Total = dades_lasioglossum$Abundancia_Total)
rownames(dades_lasioglossum_ok) <- rownames(dades_lasioglossum)

# View(dades_lasioglossum_ok)

# Extreure columnes de zona i tractament
lasioglossums_zona_tractament <- dades_lasioglossum_ok %>%
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  mutate(Tractament = str_extract(rownames(.), "(?<=_)[A-Z]+(?=_\\d)"))

# View(lasioglossums_zona_tractament)

Model lineal mixt (tractament com a variable explicativa, i zona com a variable aleatòria)

# install.packages("lme4")
# install.packages("Matrix")

model_lasioglossum <- lmer(Abundancia_Total ~ Tractament + (1|Zona), data = lasioglossums_zona_tractament)
summary(model_lasioglossum)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Abundancia_Total ~ Tractament + (1 | Zona)
##    Data: lasioglossums_zona_tractament
## 
## REML criterion at convergence: 642.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1440 -0.4778 -0.1408  0.1962  7.2804 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Zona     (Intercept)  1.557   1.248   
##  Residual             35.226   5.935   
## Number of obs: 101, groups:  Zona, 3
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)   
## (Intercept)    1.926      1.172  2.051   1.643  0.23908   
## TractamentS    3.954      1.226 98.906   3.224  0.00171 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.429
plot(model_lasioglossum)

Quan el tractament és NS, l’abundància de Lasioglossums és superior a 0 (Intercept: marginalment significatiu). Quan el tractament és S, s’espera que l’abundància de lasioglossums augmenti 3,954 respecte el tractament NS, i és un resultat molt significatiu (p-valor = 0.00171). Els Lasioglossums son més abundants a les zones segades (concorda amb algunes de les espècies “indicadores” de S!)

Model lineal (Tractament com a variable explicativa, sense considerar la zona)

# Convertir la variable Tractament a factor
lasioglossums_zona_tractament$Tractament <- factor(lasioglossums_zona_tractament$Tractament)

# Ajustar el model lineal
model <- lm(Abundancia_Total ~ Tractament, data = lasioglossums_zona_tractament)
summary(model)
## 
## Call:
## lm(formula = Abundancia_Total ~ Tractament, data = lasioglossums_zona_tractament)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.421 -2.619 -1.421  1.381 43.579 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.619      0.754   3.473 0.000764 ***
## TractamentS    3.802      1.229   3.093 0.002576 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.985 on 99 degrees of freedom
## Multiple R-squared:  0.08811,    Adjusted R-squared:  0.0789 
## F-statistic: 9.566 on 1 and 99 DF,  p-value: 0.002576
# Crear el gráfico
plot(Abundancia_Total ~ Tractament, data = lasioglossums_zona_tractament, pch = 16)

Diversitat de Shannon d’abelles (zona com a variable aleatòria)

# Crear un dataframe de la diversitat de shannon d'abelles (calculada abans)
taula_shannon_abelles <- data.frame(Shannon = shannondiv_abelles)
rownames(taula_shannon_abelles) <- rownames(data_shannon_abelles)

# Extreure dues noves columnes: tractament i zona
taula_shannon_abelles <- taula_shannon_abelles %>%
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2])

print(taula_shannon_abelles)
##                          Shannon       Zona Tractament
## Campus-juliol-2020_NS   2.611229     Campus         NS
## Campus-juliol-2020_S    1.905899     Campus          S
## Campus-Juliol-2021_NS   2.090031     Campus         NS
## Campus-Juliol-2021_S    1.806531     Campus          S
## Campus-Jun-2021_NS      2.701045     Campus         NS
## Campus-Jun-2021_S       2.080060     Campus          S
## Campus-maig-2020_NS     2.199104     Campus         NS
## Franqueses-Maig-2021_NS 2.239746 Franqueses         NS
## Franqueses-Maig-2021_S  1.560710 Franqueses          S
## Sabadell-Juny-2021_NS   2.045357   Sabadell         NS
## Sabadell-Juny-2021_S    1.884985   Sabadell          S
# Fer el model
model_shannonabelles <- lmer(Shannon ~ Tractament + (1|Zona), data = taula_shannon_abelles)
summary(model_shannonabelles)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon ~ Tractament + (1 | Zona)
##    Data: taula_shannon_abelles
## 
## REML criterion at convergence: 3.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.12461 -0.79509 -0.01592  0.61598  1.56288 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Zona     (Intercept) 0.009103 0.09541 
##  Residual             0.051690 0.22736 
## Number of obs: 11, groups:  Zona, 3
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   2.2831     0.1121  3.4370  20.373 0.000105 ***
## TractamentS  -0.4605     0.1378  7.5671  -3.342 0.011055 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.570

Quan el tractament és NS, el valor de l’índex de Shannon és significativament superior a 0 (Intercept). Quan el tractament és S, s’espera que l’índex de Shannon disminueixi 0,4605 respecte el tractament NS. És un resultat significatiu (p-valor = 0,011055). Per tant, la diversitat de Shannon en abelles és significativament inferior en el tractament segat.

# Plot del model
plot(fitted(model_shannonabelles), resid(model_shannonabelles), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Riquesa total

# Preparar les dades
riquesa_aux <- as.data.frame(riquesa)
rownames(riquesa_aux) <- paste(riquesa_aux$Mostreig, riquesa_aux$Tractament, sep = "_")
riquesa_aux <- subset(riquesa_aux, select = -Mostreig)

# Extreure columna de la zona
taula_riquesa <- riquesa_aux %>%
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1])

print(taula_riquesa)
##                         Tractament riquesa       Zona
## Campus-Juliol-2021_NS           NS      20     Campus
## Campus-Juliol-2021_S             S      21     Campus
## Campus-Jun-2021_NS              NS      45     Campus
## Campus-Jun-2021_S                S      34     Campus
## Campus-juliol-2020_NS           NS      39     Campus
## Campus-juliol-2020_S             S      24     Campus
## Campus-maig-2020_NS             NS      45     Campus
## Franqueses-Maig-2021_NS         NS      20 Franqueses
## Franqueses-Maig-2021_S           S      11 Franqueses
## Sabadell-Juny-2021_NS           NS      16   Sabadell
## Sabadell-Juny-2021_S             S      12   Sabadell
# Fer el model
model_riquesa <- lmer(riquesa ~ Tractament + (1|Zona), data = taula_riquesa)
summary(model_riquesa)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: riquesa ~ Tractament + (1 | Zona)
##    Data: taula_riquesa
## 
## REML criterion at convergence: 70.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7954 -0.4064 -0.1774  0.6667  1.1012 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Zona     (Intercept) 86.50    9.301   
##  Residual             74.49    8.631   
## Number of obs: 11, groups:  Zona, 3
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)  
## (Intercept)   26.171      6.627  2.885   3.949   0.0311 *
## TractamentS   -9.501      5.237  7.267  -1.814   0.1109  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.373

TractamentS: quan el tractament és S la riquesa disminueix (-9.501), però no és un resultat significatiu.

# Plot del model
plot(fitted(model_riquesa), resid(model_riquesa), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Riquesa PER TRAMPA d’abelles al Campus (data de mostreig com a factor aleatori)

# Preparar les dades
dades_abellescampus <- dades_abelles[dades_abelles$Zona == "Campus", ]

riquesatrampa_abellescampus <- dades_abellescampus %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Riquesa = n_distinct(ID),.groups = "drop")

riquesa_abellescampus <- as.data.frame(riquesatrampa_abellescampus)
rownames(riquesa_abellescampus) <- paste(riquesa_abellescampus$Mostreig, riquesa_abellescampus$Tractament, riquesa_abellescampus$Trampa, sep = "_")
riquesa_abellescampus <- subset(riquesa_abellescampus, select = -Mostreig)

# Extreure columna de la data de mostreig:
riquesa_abellescampus <- riquesa_abellescampus %>%
  mutate(Data_Mostreig = str_extract(rownames(.), "(?<=-)[^-_]+-[0-9]+"))

head(riquesa_abellescampus)
##                          Tractament Trampa Riquesa Data_Mostreig
## Campus-Juliol-2021_NS_2          NS      2       1   Juliol-2021
## Campus-Juliol-2021_NS_4          NS      4       1   Juliol-2021
## Campus-Juliol-2021_NS_6          NS      6       2   Juliol-2021
## Campus-Juliol-2021_NS_8          NS      8       3   Juliol-2021
## Campus-Juliol-2021_NS_10         NS     10       1   Juliol-2021
## Campus-Juliol-2021_NS_12         NS     12       3   Juliol-2021
# Fer el model
model_riquesa_abellescampus <- lmer(Riquesa ~ Tractament + (1|Data_Mostreig), data = riquesa_abellescampus)
summary(model_riquesa_abellescampus)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Riquesa ~ Tractament + (1 | Data_Mostreig)
##    Data: riquesa_abellescampus
## 
## REML criterion at convergence: 265
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8026 -0.7445 -0.1710  0.4959  2.9434 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  Data_Mostreig (Intercept) 0.6885   0.8297  
##  Residual                  4.4396   2.1070  
## Number of obs: 61, groups:  Data_Mostreig, 3
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)   
## (Intercept)   4.0378     0.6041  2.8500   6.684  0.00801 **
## TractamentS  -0.3864     0.5427 57.1316  -0.712  0.47932   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.414

En el tractament S la riquesa d’abelles per trampa disminueix, però no és significatiu.

# Plot del model
plot(fitted(model_riquesa_abellescampus), resid(model_riquesa_abellescampus), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Riquesa PER MOSTREIG d’abelles al Campus (data de mostreig com a factor aleatori)

# Preparar les dades
riquesa_aux <- dades_abellescampus %>%
  group_by(Mostreig, Tractament) %>%
  summarise(riquesa = n_distinct(ID))
## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.
riquesa_aux <- as.data.frame(riquesa_aux)
rownames(riquesa_aux) <- paste(riquesa_aux$Mostreig, riquesa_aux$Tractament, sep = "_")
riquesa_aux <- subset(riquesa_aux, select = -Mostreig)

# Extreure columna de la data de mostreig:
riquesa_aux <- riquesa_aux %>%
  mutate(Data_Mostreig = str_extract(rownames(.), "(?<=-)[^-_]+-[0-9]+"))

riquesa_aux
##                       Tractament riquesa Data_Mostreig
## Campus-Juliol-2021_NS         NS      12   Juliol-2021
## Campus-Juliol-2021_S           S      17   Juliol-2021
## Campus-Jun-2021_NS            NS      25      Jun-2021
## Campus-Jun-2021_S              S      19      Jun-2021
## Campus-juliol-2020_NS         NS      27   juliol-2020
## Campus-juliol-2020_S           S      14   juliol-2020
# Fer el model
model_riquesa_abellescampus2 <- lmer(riquesa ~ Tractament + (1|Data_Mostreig), data = riquesa_aux)
## boundary (singular) fit: see help('isSingular')
summary(model_riquesa_abellescampus2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: riquesa ~ Tractament + (1 | Data_Mostreig)
##    Data: riquesa_aux
## 
## REML criterion at convergence: 27.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5484 -0.3180  0.2212  0.5530  0.9401 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  Data_Mostreig (Intercept)  0.00    0.000   
##  Residual                  36.33    6.028   
## Number of obs: 6, groups:  Data_Mostreig, 3
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)   
## (Intercept)   21.333      3.480  4.000   6.130  0.00359 **
## TractamentS   -4.667      4.922  4.000  -0.948  0.39672   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.707
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

TractamentS: quan el tractament és S la riquesa disminueix (en 5.583), però NO és un resultat significatiu (p-valor = 0.24677).

# Plot del model
plot(fitted(model_riquesa_abellescampus2), resid(model_riquesa_abellescampus2), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Diversitat (Shannon) d’abelles al Campus (data de mostreig com a factor aleatori)

# Crear taula auxiliar: abundàncies d'abelles per ID, per cada mostreig del Campus

# Filtrar dades per Campus
dades_abellescampus <- dades_finals[dades_finals$Zona == "Campus", ]

# Crear taula de dades d'abundància només per abelles
abundanciabelles_ID2 <- dades_abellescampus %>% 
  mutate(Mostreig = paste(Zona, Mes, Any, sep = "-")) %>% 
  filter(Ordre == 'HymenopteraAbella') %>% 
  group_by(Mostreig, Tractament, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Mostreig', 'Tractament'. You can override
## using the `.groups` argument.
# View(abundanciabelles_ID)

#Crear una taula auxiliar de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux2 <- cast(abundanciabelles_ID2, Mostreig + Tractament ~ ID, value='Numero_total', FUN=mean)
taula_aux2 <- as.data.frame(taula_aux2)
taula_aux2[is.na(taula_aux2)] <- 0

#Agrupar Mostreig i tractament en una nova variable (Agrupacio)
data_shannon_abellescampus <- taula_aux2 %>% 
  mutate(Agrupacio = paste(Mostreig, Tractament, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament)

#Definir els rownames amb la variable Agrupacio (mostreig i tractament)
rownames(data_shannon_abellescampus) <- data_shannon_abellescampus$Agrupacio
data_shannon_abellescampus <- data_shannon_abellescampus %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)

# head(data_shannon_abellescampus)
# Càlcul de l'índex de Shannon i preparació del dataframe
shannondiv_abellescampus <- diversity(data_shannon_abellescampus)

# Crear un dataframe de la diversitat de shannon d'abelles
taula_shannon_abellescampus <- data.frame(Shannon = shannondiv_abellescampus)
rownames(taula_shannon_abellescampus) <- rownames(data_shannon_abellescampus)

# Extreure dues noves columnes: tractament i data mostreig
taula_shannon_abellescampus <- taula_shannon_abellescampus %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% 
  mutate(Data_Mostreig = str_extract(rownames(.), "(?<=-)[^-_]+-[0-9]+"))

print(taula_shannon_abellescampus)
##                        Shannon Tractament Data_Mostreig
## Campus-juliol-2020_NS 2.611229         NS   juliol-2020
## Campus-juliol-2020_S  1.905899          S   juliol-2020
## Campus-Juliol-2021_NS 2.090031         NS   Juliol-2021
## Campus-Juliol-2021_S  1.806531          S   Juliol-2021
## Campus-Jun-2021_NS    2.701045         NS      Jun-2021
## Campus-Jun-2021_S     2.080060          S      Jun-2021
## Campus-maig-2020_NS   2.199104         NS     maig-2020
# Fer el model
model_shannonabellescampus <- lmer(Shannon ~ Tractament + (1|Data_Mostreig), data = taula_shannon_abellescampus)
summary(model_shannonabellescampus)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon ~ Tractament + (1 | Data_Mostreig)
##    Data: taula_shannon_abellescampus
## 
## REML criterion at convergence: 1.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.03570 -0.47827  0.03651  0.57598  0.80378 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  Data_Mostreig (Intercept) 0.03686  0.1920  
##  Residual                  0.02461  0.1569  
## Number of obs: 7, groups:  Data_Mostreig, 4
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   2.4004     0.1240  3.9717  19.363 4.42e-05 ***
## TractamentS  -0.5097     0.1248  2.2656  -4.083   0.0444 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.398

TractamentS: quan el tractament és S, s’espera que l’índex de Shannon disminueixi 0,5097 respecte el tractament NS. És un resultat significatiu! (p-valor= 0.0444).

# Plot del model
plot(fitted(model_shannonabellescampus), resid(model_shannonabellescampus), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

MODELS LINEALS D’ANÀLISIS FUNCIONALS

Inter-tegular span (ITS) (model lineal mixt)

Abelles totals, amb Zona com a factor aleatori

#Obrir la base de dades amb la columna de les mesures de la ITS i filtrar dades per abelles
abelles_ITS <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]

head(abelles_ITS)
## # A tibble: 6 × 26
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6       159 Campus Juliol         7  2021     10 NS         <NA>         
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Model lineal mixt ITS (mida)
abelles_ITS_filtrades <- abelles_ITS[!is.na(abelles_ITS$Inter_tegular_span),]
model_abelles_ITS <- lmer(`Inter_tegular_span` ~ Tractament + (1|Zona), data = abelles_ITS_filtrades)
summary(model_abelles_ITS)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Inter_tegular_span ~ Tractament + (1 | Zona)
##    Data: abelles_ITS_filtrades
## 
## REML criterion at convergence: 1453.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3928 -0.4191 -0.2937  0.1489  5.4311 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Zona     (Intercept) 0.03073  0.1753  
##  Residual             0.48378  0.6955  
## Number of obs: 683, groups:  Zona, 3
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   1.78796    0.11569   2.44978  15.455  0.00167 ** 
## TractamentS  -0.51521    0.05431 676.60529  -9.487  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.217

La ITS de les abelles disminueix -0.51521 quan el tractament és S (abelles més petites quan està segat). És significatiu (p-valor = 2e-16)

# Plot del model
# Prediccions del model sobre les dades:
predicted_values <- fitted(model_abelles_ITS)
observed_values <- abelles_ITS[!is.na(abelles_ITS$Inter_tegular_span),]$`Inter_tegular_span`

# Crear el data frame pel gràfic
plot_data <- data.frame(Observed = observed_values, Predicted = predicted_values)

# Graficar
ggplot(plot_data, aes(x = Predicted, y = Observed)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, color = "green", linetype = "dashed") +  # Línea de referencia de 45 grados
  labs(x = "Valors predits", y = "Valors observats", title = "Valors observats vs. predits") +
  theme_minimal()

Visualitzar gràficament amb boxplot

ggplot(abelles_ITS_filtrades, aes(Tractament, Inter_tegular_span, fill = Tractament)) +
  geom_boxplot() +
  facet_wrap(~Zona) +
  labs(x = "Tratament", y = "ITS") +
  scale_fill_manual(values = c("S" = "#ADD8E6", "NS" = "#FFC0CB"))

Abelles del Campus, amb data_mostreig com a factor aleatori

# Filtrar les dades només per abelles del Campus
abellescampus_ITS <- abelles_ITS[abelles_ITS$Zona == "Campus",]

#Definir nova variable: "data_mostreig", combinació de mes i any de mostreig
abellescampus_ITS <- abellescampus_ITS %>%
  mutate(data_mostreig = paste(Mes, Any, sep = "-"))

head(abellescampus_ITS)
## # A tibble: 6 × 27
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6       159 Campus Juliol         7  2021     10 NS         <NA>         
## # ℹ 19 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>, data_mostreig <chr>
# Fer el model
model_abellescampus_ITS <- lmer(`Inter_tegular_span` ~ Tractament + (1|data_mostreig), data = abellescampus_ITS)
summary(model_abellescampus_ITS)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Inter_tegular_span ~ Tractament + (1 | data_mostreig)
##    Data: abellescampus_ITS
## 
## REML criterion at convergence: 1196.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3101 -0.5834 -0.2132  0.1382  5.4645 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  data_mostreig (Intercept) 0.04629  0.2151  
##  Residual                  0.42778  0.6540  
## Number of obs: 594, groups:  data_mostreig, 4
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   1.54523    0.11430   3.29968  13.519 0.000529 ***
## TractamentS  -0.36675    0.06534 537.59881  -5.613 3.19e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.243

Molt significatiu: la ITS de les abelles (mida) disminueix -0.36675 quan el tractament és S. Les abelles son més petites quan està segat.

# Fer el plot
# Prediccions del model sobre les dades:
predicted_values1 <- fitted(model_abellescampus_ITS)
observed_values1 <- abellescampus_ITS[!is.na(abellescampus_ITS$Inter_tegular_span),]$`Inter_tegular_span`

# Crear el data frame pel gràfic
plot_data <- data.frame(Observed = observed_values1, Predicted = predicted_values1)

# Graficar
ggplot(plot_data, aes(x = Predicted, y = Observed)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, color = "green", linetype = "dashed") +  
  labs(x = "Valors predits", y = "Valors observats", title = "Valors observats vs. predits") +
  theme_minimal()

### ATENCIÓ!!! NO ESTIC SEGURA DE SI AQUESTA REPRESENTACIÓ GRÀFICA DEL MODEL ÉS ACURADA. SON VALORS PREDITS I OBSERVATS, O SON ELS RESIDUS AJUSTATS I NOSEQ???? VEURE GRÀFICS DELS ALTRES MODELS ANTERIORS.

ITS abelles per espècies totals:

# Crear nova base de dades: només una espècie per cada tractament.
especies_ITS <- abelles_ITS %>%
  distinct(ID, Tractament, .keep_all = TRUE)
# Dialectus -- dialectus
head(especies_ITS)
## # A tibble: 6 × 26
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6       159 Campus Juliol         7  2021     10 NS         <NA>         
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Fer el model (lineal, no mixt: no podem considerar la zona com a factor aleatori ja que hem eliminat espècimens!)
model_especies_ITS <- lm(`Inter_tegular_span` ~ Tractament, data = especies_ITS)
summary(model_especies_ITS)
## 
## Call:
## lm(formula = Inter_tegular_span ~ Tractament, data = especies_ITS)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.3679 -0.7225 -0.3465  0.6686  3.5533 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.1199     0.1603  13.225   <2e-16 ***
## TractamentS  -0.2894     0.2581  -1.121    0.266    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.051 on 68 degrees of freedom
##   (35 observations deleted due to missingness)
## Multiple R-squared:  0.01815,    Adjusted R-squared:  0.003711 
## F-statistic: 1.257 on 1 and 68 DF,  p-value: 0.2662

No significatiu

ITS abelles per espècies del campus

# Crear nova base de dades: només una espècie per cada tractament.
especiescampus_ITS <- abellescampus_ITS %>%
  distinct(ID, Tractament, .keep_all = TRUE)
# Dialectus -- dialectus
View(especiescampus_ITS)
# Fer el model (lineal, no mixt: no podem considerar la zona com a factor aleatori ja que hem eliminat espècimens!)
model_especiescampus_ITS <- lm(`Inter_tegular_span` ~ Tractament, data = especiescampus_ITS)
summary(model_especiescampus_ITS)
## 
## Call:
## lm(formula = Inter_tegular_span ~ Tractament, data = especiescampus_ITS)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2394 -0.6734 -0.3183  0.6081  2.8422 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.9914     0.1667  11.949   <2e-16 ***
## TractamentS  -0.2301     0.2582  -0.891    0.377    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9859 on 58 degrees of freedom
##   (33 observations deleted due to missingness)
## Multiple R-squared:  0.01351,    Adjusted R-squared:  -0.003501 
## F-statistic: 0.7941 on 1 and 58 DF,  p-value: 0.3765

No significatiu.

Substrat de nidificació (GLMM binomial)

Nest_type en funció del tractament de les abelles totals (zona = factor aleatori)

abelles_ITS<-abelles_ITS[!is.na(abelles_ITS$Nesting_type),]

abelles_ITS["Nius"] = 0
abelles_ITS$Nius[abelles_ITS["Nesting_type"]=="S"]="S"
abelles_ITS$Nius[abelles_ITS["Nesting_type"]!="S"]="Not_S"

# Taula de contingència
taula_contingencia1 <- table(abelles_ITS$Tractament, abelles_ITS$Nius)
print(taula_contingencia1)
##     
##      Not_S   S
##   NS    71 324
##   S     14 291
# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia1)
ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
  labs(x = "Tractament", y = "Nº", fill = "Nesting_type") +
  ggtitle("Substrats de nidificació per tractament")

# Model:
abelles_ITS$Nius <- factor(abelles_ITS$Nius)
# Nomes Zona o tambe mostreig, com a factors aleatoris
model_abelles_nest <- glmer(Nius ~ Tractament + (1|Zona) + (1|Mostreig), data = abelles_ITS, family = binomial)
## boundary (singular) fit: see help('isSingular')
summary(model_abelles_nest)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Nius ~ Tractament + (1 | Zona) + (1 | Mostreig)
##    Data: abelles_ITS
## 
##      AIC      BIC   logLik deviance df.resid 
##    493.3    511.5   -242.6    485.3      696 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7820  0.2091  0.2282  0.4715  0.5146 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  Mostreig (Intercept) 3.795e-02 0.1947978
##  Zona     (Intercept) 1.277e-10 0.0000113
## Number of obs: 700, groups:  Mostreig, 6; Zona, 3
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.4843     0.1645   9.024  < 2e-16 ***
## TractamentS   1.6262     0.3450   4.714 2.43e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.442
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

Significatiu! Les abelles del tractament S tenen més probabilitat (1.6262) de fer nius al sòl que de fer nius a altres llocs que no siguin el sòl. És un resultat significatiu! (p-valor = 2.43e-06).

# Plot del model
plot(fitted(model_abelles_nest), resid(model_abelles_nest), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Nest_type en funció del tractament de les abelles del Campus (data_mostreig = factor aleatori)

abellescampus_ITS<-abellescampus_ITS[!is.na(abellescampus_ITS$Nesting_type),]

abellescampus_ITS["Nius"] = 0
abellescampus_ITS$Nius[abellescampus_ITS["Nesting_type"]=="S"]="S"
abellescampus_ITS$Nius[abellescampus_ITS["Nesting_type"]!="S"]="Not_S"

# Taula de contingència
taula_contingencia2 <- table(abellescampus_ITS$Tractament, abellescampus_ITS$Nius)
print(taula_contingencia2)
##     
##      Not_S   S
##   NS    63 291
##   S     13 243
# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia2)

ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
  labs(x = "Tractament", y = "Nº", fill = "Nesting_type") +
  ggtitle("Substrats de nidificació per tractament")

# Model:
abellescampus_ITS$Nesting_type <- factor(abellescampus_ITS$Nius)
model_abellescampus_nest <- glmer(Nesting_type ~ Tractament + (1|data_mostreig), data = abellescampus_ITS, family = binomial)

summary(model_abellescampus_nest)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Nesting_type ~ Tractament + (1 | data_mostreig)
##    Data: abellescampus_ITS
## 
##      AIC      BIC   logLik deviance df.resid 
##    440.2    453.4   -217.1    434.2      607 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5162  0.2291  0.2360  0.4689  0.4997 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  data_mostreig (Intercept) 0.02712  0.1647  
## Number of obs: 610, groups:  data_mostreig, 4
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.4968     0.1764   8.484  < 2e-16 ***
## TractamentS   1.5004     0.3860   3.887 0.000101 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.512

Significatiu! Les abelles del tractament S tenen més probabilitat (1.5004) de fer nius al sòl que de fer nius a altres llocs que no siguin el sòl. És un resultat significatiu! (p-valor = 0.000101).

# Plot del model
plot(fitted(model_abellescampus_nest), resid(model_abellescampus_nest), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Lèctia (GLMM binomial)

Lecty en funció del tractament de les abelles totals (Zona = factor aleatori)

# Taula de contingència
taula_contingencia3 <- table(abelles_ITS$Tractament, abelles_ITS$Lecty)
print(taula_contingencia3)
##     
##      Oligolectic Polylectic
##   NS         142        249
##   S           38        263
# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia3)

ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
  labs(x = "Tractament", y = "Nº", fill = "Lecty") +
  ggtitle("Lecty per tractament")

# Model:
abelles_ITS$Lecty <- factor(abelles_ITS$Lecty)
model_abelles_lecty <- glmer(Lecty ~ Tractament + (1|Zona), data = abelles_ITS, family = binomial)

summary(model_abelles_lecty)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Lecty ~ Tractament + (1 | Zona)
##    Data: abelles_ITS
## 
##      AIC      BIC   logLik deviance df.resid 
##    730.0    743.6   -362.0    724.0      689 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0678 -0.6616  0.3260  0.7371  1.5115 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Zona   (Intercept) 0.7457   0.8635  
## Number of obs: 692, groups:  Zona, 3
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.2956     0.5350   0.553    0.581    
## TractamentS   1.6319     0.2235   7.302 2.84e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.132

Les abelles tenen una major probabilitat de ser polilèctiques quan estan sotmeses al tractament S. És un resultat molt significatiu (p-valor=), la qual cosa suggereix un fort efecte del tractament en la lecticitat de les abelles.

# Plot del model
plot(fitted(model_abelles_lecty), resid(model_abelles_lecty), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Lecty en funció del tractament de les abelles del campus (data_mostreig = factor aleatori)

# Taula de contingència
taula_contingencia4 <- table(abellescampus_ITS$Tractament, abellescampus_ITS$Lecty)
print(taula_contingencia4)
##     
##      Oligolectic Polylectic
##   NS         127        223
##   S           20        232
# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia4)

ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
  labs(x = "Tractament", y = "Nº", fill = "Lecty") +
  ggtitle("Lecty per tractament")

# Model:
abellescampus_ITS$Lecty <- factor(abellescampus_ITS$Lecty)
model_abellescampus_lecty <- glmer(Lecty ~ Tractament + (1|data_mostreig), data = abellescampus_ITS, family = binomial)

summary(model_abellescampus_lecty)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Lecty ~ Tractament + (1 | data_mostreig)
##    Data: abellescampus_ITS
## 
##      AIC      BIC   logLik deviance df.resid 
##    468.8    482.0   -231.4    462.8      599 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2618  0.0518  0.1706  0.5796  1.2245 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  data_mostreig (Intercept) 6.291    2.508   
## Number of obs: 602, groups:  data_mostreig, 4
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   2.2170     1.3182   1.682   0.0926 .
## TractamentS   0.6375     0.3271   1.949   0.0513 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.091

Les abelles tenen una major probabilitat de ser polilèctiques quan estan sotmeses al tractament S. És un resultat marginalment significatiu (p-valor = 0.0513), la qual cosa suggereix un cert efecte del tractament en la lecticitat de les abelles.

# Plot del model
plot(fitted(model_abellescampus_lecty), resid(model_abellescampus_lecty), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")